Add newtonian line search

.pre-commit-config.yaml

@@ -10,8 +10,8 @@ repos:
       - id: check-yaml
         exclude: ^data/
 
-  - repo: https://github.com/pre-commit/mirrors-clang-format
-    rev: v18.1.8
+  - repo: https://github.com/pocc/pre-commit-hooks
+    rev: v1.3.5
     hooks:
       - id: clang-format
         name: Fix *.c,*.cc,*.cpp,*.h,*.hpp

python/csrc/kintera.cpp

@@ -28,13 +28,13 @@ extern std::vector<double> species_sref_R;
 
 }  // namespace kintera
 
-void bind_thermo(py::module &m);
-void bind_constants(py::module &m);
-void bind_kinetics(py::module &m);
-void bind_photolysis(py::module &m);
-void bind_diffusion(py::module &m);
-void bind_sparse_solver(py::module &m);
-void bind_equilibrium(py::module &m);
+void bind_thermo(py::module& m);
+void bind_constants(py::module& m);
+void bind_kinetics(py::module& m);
+void bind_photolysis(py::module& m);
+void bind_diffusion(py::module& m);
+void bind_sparse_solver(py::module& m);
+void bind_equilibrium(py::module& m);
 
 PYBIND11_MODULE(kintera, m) {
   m.attr("__name__") = "kintera";
@@ -55,7 +55,7 @@ PYBIND11_MODULE(kintera, m) {
 
   pySpeciesThermo.def(py::init<>())
       .def("__repr__",
-           [](const kintera::SpeciesThermo &self) {
+           [](const kintera::SpeciesThermo& self) {
              return fmt::format("SpeciesThermo({})", self);
            })
       .def("species", &kintera::SpeciesThermoImpl::species)
@@ -70,9 +70,9 @@ PYBIND11_MODULE(kintera, m) {
   auto pyReaction = py::class_<kintera::Reaction>(m, "Reaction");
 
   pyReaction.def(py::init<>())
-      .def(py::init<const std::string &>())
+      .def(py::init<const std::string&>())
       .def("__repr__",
-           [](const kintera::Reaction &self) {
+           [](const kintera::Reaction& self) {
              return fmt::format("Reaction({})", self);
            })
       .def("equation", &kintera::Reaction::equation)
@@ -220,7 +220,7 @@ PYBIND11_MODULE(kintera, m) {
                      &kintera::KBTitanReactionReport::unsupported_reaction_ids);
 
   m.def("classify_kinetics_base_titan_reactions",
-        py::overload_cast<std::string const &, std::string const &>(
+        py::overload_cast<std::string const&, std::string const&>(
             &kintera::classify_kinetics_base_titan_reactions),
         py::arg("pun_path"), py::arg("run_input_path"));
 
@@ -232,47 +232,44 @@ PYBIND11_MODULE(kintera, m) {
   bind_sparse_solver(m);
   bind_equilibrium(m);
 
-  m.def("species_names", []() -> const std::vector<std::string> & {
+  m.def("species_names", []() -> const std::vector<std::string>& {
     return kintera::species_names;
   });
 
-  m.def("set_species_names", [](const std::vector<std::string> &names) {
+  m.def("set_species_names", [](const std::vector<std::string>& names) {
     kintera::species_names = names;
     return kintera::species_names;
   });
 
-  m.def("species_weights", []() -> const std::vector<double> & {
+  m.def("species_weights", []() -> const std::vector<double>& {
     return kintera::species_weights;
   });
 
-  m.def("set_species_weights", [](const std::vector<double> &weights) {
+  m.def("set_species_weights", [](const std::vector<double>& weights) {
     kintera::species_weights = weights;
     return kintera::species_weights;
   });
 
-  m.def("species_cref_R", []() -> const std::vector<double> & {
-    return kintera::species_cref_R;
-  });
+  m.def("species_cref_R",
+        []() -> const std::vector<double>& { return kintera::species_cref_R; });
 
-  m.def("set_species_cref_R", [](const std::vector<double> &cref_R) {
+  m.def("set_species_cref_R", [](const std::vector<double>& cref_R) {
     kintera::species_cref_R = cref_R;
     return kintera::species_cref_R;
   });
 
-  m.def("species_uref_R", []() -> const std::vector<double> & {
-    return kintera::species_uref_R;
-  });
+  m.def("species_uref_R",
+        []() -> const std::vector<double>& { return kintera::species_uref_R; });
 
-  m.def("set_species_uref_R", [](const std::vector<double> &uref_R) {
+  m.def("set_species_uref_R", [](const std::vector<double>& uref_R) {
     kintera::species_uref_R = uref_R;
     return kintera::species_uref_R;
   });
 
-  m.def("species_sref_R", []() -> const std::vector<double> & {
-    return kintera::species_sref_R;
-  });
+  m.def("species_sref_R",
+        []() -> const std::vector<double>& { return kintera::species_sref_R; });
 
-  m.def("set_species_sref_R", [](const std::vector<double> &sref_R) {
+  m.def("set_species_sref_R", [](const std::vector<double>& sref_R) {
     kintera::species_sref_R = sref_R;
     return kintera::species_sref_R;
   });

python/csrc/pyequilibrium.cpp

@@ -9,14 +9,14 @@
 
 namespace py = pybind11;
 
-void bind_equilibrium(py::module &m) {
+void bind_equilibrium(py::module& m) {
   auto pyOptions =
       py::class_<kintera::EquilibriumOptionsImpl, kintera::EquilibriumOptions>(
           m, "EquilibriumOptions");
 
   pyOptions.def(py::init<>())
       .def("__repr__",
-           [](kintera::EquilibriumOptions const &self) {
+           [](kintera::EquilibriumOptions const& self) {
              std::stringstream ss;
              self->report(ss);
              return ss.str();

src/equilibrium/equilibrium.cpp

@@ -11,7 +11,7 @@
 
 namespace kintera {
 
-void EquilibriumOptionsImpl::report(std::ostream &os) const {
+void EquilibriumOptionsImpl::report(std::ostream& os) const {
   os << "-- equilibrium options --\n"
      << "* components = " << components().size() << "\n"
      << "* phases = " << phases().size() << "\n"
@@ -42,15 +42,15 @@ void EquilibriumOptionsImpl::validate() const {
   TORCH_CHECK(ftol() > 0., "ftol must be positive");
   TORCH_CHECK(mole_floor() >= 0., "mole_floor must be nonnegative");
 
-  for (auto const &equation : reactions()) {
+  for (auto const& equation : reactions()) {
     Reaction reaction(equation);
-    for (auto const &[name, coefficient] : reaction.reactants()) {
+    for (auto const& [name, coefficient] : reaction.reactants()) {
       (void)coefficient;
       TORCH_CHECK(std::find(components().begin(), components().end(), name) !=
                       components().end(),
                   "reaction references unknown component: ", name);
     }
-    for (auto const &[name, coefficient] : reaction.products()) {
+    for (auto const& [name, coefficient] : reaction.products()) {
       (void)coefficient;
       TORCH_CHECK(std::find(components().begin(), components().end(), name) !=
                       components().end(),
@@ -59,7 +59,7 @@ void EquilibriumOptionsImpl::validate() const {
   }
 }
 
-EquilibriumTPImpl::EquilibriumTPImpl(EquilibriumOptions const &options_)
+EquilibriumTPImpl::EquilibriumTPImpl(EquilibriumOptions const& options_)
     : options(options_) {
   reset();
 }
@@ -68,7 +68,7 @@ void EquilibriumTPImpl::reset() {
   options->validate();
   std::vector<Reaction> reactions;
   reactions.reserve(options->reactions().size());
-  for (auto const &equation : options->reactions())
+  for (auto const& equation : options->reactions())
     reactions.emplace_back(equation);
   stoich = register_buffer(
       "stoich", generate_stoichiometry_matrix(reactions, options->components())
@@ -82,7 +82,7 @@ void EquilibriumTPImpl::reset() {
                               torch::tensor(options->phase_ids(), torch::kInt));
 }
 
-void EquilibriumTPImpl::pretty_print(std::ostream &os) const {
+void EquilibriumTPImpl::pretty_print(std::ostream& os) const {
   os << "EquilibriumTP(components=" << options->components().size()
      << ", phases=" << options->phases().size() << ")";
 }

src/equilibrium/equilibrium.hpp

@@ -18,9 +18,9 @@ struct EquilibriumOptionsImpl final {
   }
 
   static std::shared_ptr<EquilibriumOptionsImpl> from_yaml(
-      std::string const &filename, bool verbose = false);
+      std::string const& filename, bool verbose = false);
 
-  void report(std::ostream &os) const;
+  void report(std::ostream& os) const;
   void validate() const;
 
   ADD_ARG(std::vector<std::string>, components);
@@ -43,10 +43,10 @@ class EquilibriumTPImpl : public torch::nn::Cloneable<EquilibriumTPImpl> {
   torch::Tensor phase_ids;
 
   EquilibriumTPImpl() : options(EquilibriumOptionsImpl::create()) {}
-  explicit EquilibriumTPImpl(EquilibriumOptions const &options_);
+  explicit EquilibriumTPImpl(EquilibriumOptions const& options_);
 
   void reset() override;
-  void pretty_print(std::ostream &os) const override;
+  void pretty_print(std::ostream& os) const override;
 
   //! Solve fixed-temperature, fixed-pressure multiphase equilibrium.
   /*!

src/equilibrium/equilibrium_dispatch.cpp

@@ -8,8 +8,8 @@
 
 namespace kintera {
 
-void call_equilibrium_cpu(at::TensorIterator &iter, at::Tensor const &stoich,
-                          at::Tensor const &phase_ids, int nphase,
+void call_equilibrium_cpu(at::TensorIterator& iter, at::Tensor const& stoich,
+                          at::Tensor const& phase_ids, int nphase,
                           int gas_phase, double standard_pressure, double ftol,
                           double mole_floor, int max_iter) {
   int grain_size = std::max<int64_t>(1, iter.numel() / at::get_num_threads());
@@ -19,15 +19,15 @@ void call_equilibrium_cpu(at::TensorIterator &iter, at::Tensor const &stoich,
     auto stoich_ptr = stoich.data_ptr<scalar_t>();
     auto phase_ptr = phase_ids.data_ptr<int>();
     iter.for_each(
-        [&](char **data, const int64_t *strides, int64_t n) {
+        [&](char** data, const int64_t* strides, int64_t n) {
           for (int64_t i = 0; i < n; ++i) {
-            auto gain = reinterpret_cast<scalar_t *>(data[0] + i * strides[0]);
-            auto diag = reinterpret_cast<scalar_t *>(data[1] + i * strides[1]);
-            auto out = reinterpret_cast<scalar_t *>(data[2] + i * strides[2]);
-            auto temp = reinterpret_cast<scalar_t *>(data[3] + i * strides[3]);
-            auto pres = reinterpret_cast<scalar_t *>(data[4] + i * strides[4]);
-            auto moles = reinterpret_cast<scalar_t *>(data[5] + i * strides[5]);
-            auto log_k = reinterpret_cast<scalar_t *>(data[6] + i * strides[6]);
+            auto gain = reinterpret_cast<scalar_t*>(data[0] + i * strides[0]);
+            auto diag = reinterpret_cast<scalar_t*>(data[1] + i * strides[1]);
+            auto out = reinterpret_cast<scalar_t*>(data[2] + i * strides[2]);
+            auto temp = reinterpret_cast<scalar_t*>(data[3] + i * strides[3]);
+            auto pres = reinterpret_cast<scalar_t*>(data[4] + i * strides[4]);
+            auto moles = reinterpret_cast<scalar_t*>(data[5] + i * strides[5]);
+            auto log_k = reinterpret_cast<scalar_t*>(data[6] + i * strides[6]);
             phase_equilibrate_tp(gain, diag, out, *temp, *pres, moles, log_k,
                                  stoich_ptr, phase_ptr, nspecies, nreaction,
                                  nphase, gas_phase,

src/equilibrium/equilibrium_dispatch.hpp

@@ -5,9 +5,9 @@
 
 namespace at::native {
 
-using equilibrium_fn = void (*)(at::TensorIterator &iter,
-                                at::Tensor const &stoich,
-                                at::Tensor const &phase_ids, int nphase,
+using equilibrium_fn = void (*)(at::TensorIterator& iter,
+                                at::Tensor const& stoich,
+                                at::Tensor const& phase_ids, int nphase,
                                 int gas_phase, double standard_pressure,
                                 double ftol, double mole_floor, int max_iter);
 

src/equilibrium/equilibrium_options.cpp

@@ -13,7 +13,7 @@
 namespace kintera {
 
 EquilibriumOptions EquilibriumOptionsImpl::from_yaml(
-    std::string const &filename, bool verbose) {
+    std::string const& filename, bool verbose) {
   auto config = YAML::LoadFile(find_resource(filename));
   TORCH_CHECK(config["phases"], "equilibrium YAML requires 'phases'");
   TORCH_CHECK(config["species"], "equilibrium YAML requires 'species'");
@@ -23,7 +23,7 @@ EquilibriumOptions EquilibriumOptionsImpl::from_yaml(
   std::set<std::string> seen_components;
   int gas_phase_count = 0;
 
-  for (auto const &phase_node : config["phases"]) {
+  for (auto const& phase_node : config["phases"]) {
     TORCH_CHECK(phase_node["name"], "phase is missing 'name'");
     TORCH_CHECK(phase_node["species"], "phase is missing 'species'");
     int phase_id = static_cast<int>(options->phases().size());
@@ -40,7 +40,7 @@ EquilibriumOptions EquilibriumOptionsImpl::from_yaml(
       ++gas_phase_count;
     }
 
-    for (auto const &component_node : phase_node["species"]) {
+    for (auto const& component_node : phase_node["species"]) {
       auto component = component_node.as<std::string>();
       TORCH_CHECK(seen_components.insert(component).second,
                   "component appears in more than one phase: ", component);
@@ -53,26 +53,26 @@ EquilibriumOptions EquilibriumOptionsImpl::from_yaml(
 
   std::map<std::string, std::map<std::string, double>> compositions;
   std::set<std::string> element_set;
-  for (auto const &species_node : config["species"]) {
+  for (auto const& species_node : config["species"]) {
     TORCH_CHECK(species_node["name"], "species is missing 'name'");
     TORCH_CHECK(species_node["composition"],
                 "species is missing 'composition'");
     auto name = species_node["name"].as<std::string>();
     TORCH_CHECK(compositions.find(name) == compositions.end(),
                 "duplicate species definition: ", name);
-    for (auto const &entry : species_node["composition"]) {
+    for (auto const& entry : species_node["composition"]) {
       auto element = entry.first.as<std::string>();
       compositions[name][element] = entry.second.as<double>();
       element_set.insert(element);
     }
   }
-  for (auto const &component : options->components()) {
+  for (auto const& component : options->components()) {
     TORCH_CHECK(compositions.find(component) != compositions.end(),
                 "phase component has no species definition: ", component);
   }
 
   std::vector<Reaction> reactions;
-  for (auto const &reaction_node : config["reactions"]) {
+  for (auto const& reaction_node : config["reactions"]) {
     if (reaction_node["type"] &&
         reaction_node["type"].as<std::string>() != "equilibrium") {
       continue;
@@ -85,28 +85,28 @@ EquilibriumOptions EquilibriumOptionsImpl::from_yaml(
   TORCH_CHECK(!reactions.empty(), "equilibrium YAML contains no reactions");
 
   for (size_t j = 0; j < reactions.size(); ++j) {
-    for (auto const &[name, coefficient] : reactions[j].reactants()) {
+    for (auto const& [name, coefficient] : reactions[j].reactants()) {
       auto found = std::find(options->components().begin(),
                              options->components().end(), name);
       TORCH_CHECK(found != options->components().end(),
                   "reaction references unknown component: ", name);
       (void)coefficient;
     }
-    for (auto const &[name, coefficient] : reactions[j].products()) {
+    for (auto const& [name, coefficient] : reactions[j].products()) {
       auto found = std::find(options->components().begin(),
                              options->components().end(), name);
       TORCH_CHECK(found != options->components().end(),
                   "reaction references unknown component: ", name);
       (void)coefficient;
     }
-    for (auto const &element : element_set) {
+    for (auto const& element : element_set) {
       double balance = 0.;
-      for (auto const &[name, coefficient] : reactions[j].reactants()) {
+      for (auto const& [name, coefficient] : reactions[j].reactants()) {
         auto found = compositions.at(name).find(element);
         if (found != compositions.at(name).end())
           balance -= coefficient * found->second;
       }
-      for (auto const &[name, coefficient] : reactions[j].products()) {
+      for (auto const& [name, coefficient] : reactions[j].products()) {
         auto found = compositions.at(name).find(element);
         if (found != compositions.at(name).end())
           balance += coefficient * found->second;