diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7f8e3120..faf7fa96 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/python/csrc/kintera.cpp b/python/csrc/kintera.cpp index 0f85581f..a9d92804 100644 --- a/python/csrc/kintera.cpp +++ b/python/csrc/kintera.cpp @@ -28,13 +28,13 @@ extern std::vector 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_(m, "Reaction"); pyReaction.def(py::init<>()) - .def(py::init()) + .def(py::init()) .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( + py::overload_cast( &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 & { + m.def("species_names", []() -> const std::vector& { return kintera::species_names; }); - m.def("set_species_names", [](const std::vector &names) { + m.def("set_species_names", [](const std::vector& names) { kintera::species_names = names; return kintera::species_names; }); - m.def("species_weights", []() -> const std::vector & { + m.def("species_weights", []() -> const std::vector& { return kintera::species_weights; }); - m.def("set_species_weights", [](const std::vector &weights) { + m.def("set_species_weights", [](const std::vector& weights) { kintera::species_weights = weights; return kintera::species_weights; }); - m.def("species_cref_R", []() -> const std::vector & { - return kintera::species_cref_R; - }); + m.def("species_cref_R", + []() -> const std::vector& { return kintera::species_cref_R; }); - m.def("set_species_cref_R", [](const std::vector &cref_R) { + m.def("set_species_cref_R", [](const std::vector& cref_R) { kintera::species_cref_R = cref_R; return kintera::species_cref_R; }); - m.def("species_uref_R", []() -> const std::vector & { - return kintera::species_uref_R; - }); + m.def("species_uref_R", + []() -> const std::vector& { return kintera::species_uref_R; }); - m.def("set_species_uref_R", [](const std::vector &uref_R) { + m.def("set_species_uref_R", [](const std::vector& uref_R) { kintera::species_uref_R = uref_R; return kintera::species_uref_R; }); - m.def("species_sref_R", []() -> const std::vector & { - return kintera::species_sref_R; - }); + m.def("species_sref_R", + []() -> const std::vector& { return kintera::species_sref_R; }); - m.def("set_species_sref_R", [](const std::vector &sref_R) { + m.def("set_species_sref_R", [](const std::vector& sref_R) { kintera::species_sref_R = sref_R; return kintera::species_sref_R; }); diff --git a/python/csrc/pyequilibrium.cpp b/python/csrc/pyequilibrium.cpp index 18e53f57..fb52562a 100644 --- a/python/csrc/pyequilibrium.cpp +++ b/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_( 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(); diff --git a/src/equilibrium/equilibrium.cpp b/src/equilibrium/equilibrium.cpp index 68143636..05caefd6 100644 --- a/src/equilibrium/equilibrium.cpp +++ b/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 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() << ")"; } diff --git a/src/equilibrium/equilibrium.hpp b/src/equilibrium/equilibrium.hpp index 57200d82..12e7d1bc 100644 --- a/src/equilibrium/equilibrium.hpp +++ b/src/equilibrium/equilibrium.hpp @@ -18,9 +18,9 @@ struct EquilibriumOptionsImpl final { } static std::shared_ptr 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, components); @@ -43,10 +43,10 @@ class EquilibriumTPImpl : public torch::nn::Cloneable { 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. /*! diff --git a/src/equilibrium/equilibrium_dispatch.cpp b/src/equilibrium/equilibrium_dispatch.cpp index 167a33e1..ee4fd088 100644 --- a/src/equilibrium/equilibrium_dispatch.cpp +++ b/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(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(); auto phase_ptr = phase_ids.data_ptr(); 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(data[0] + i * strides[0]); - auto diag = reinterpret_cast(data[1] + i * strides[1]); - auto out = reinterpret_cast(data[2] + i * strides[2]); - auto temp = reinterpret_cast(data[3] + i * strides[3]); - auto pres = reinterpret_cast(data[4] + i * strides[4]); - auto moles = reinterpret_cast(data[5] + i * strides[5]); - auto log_k = reinterpret_cast(data[6] + i * strides[6]); + auto gain = reinterpret_cast(data[0] + i * strides[0]); + auto diag = reinterpret_cast(data[1] + i * strides[1]); + auto out = reinterpret_cast(data[2] + i * strides[2]); + auto temp = reinterpret_cast(data[3] + i * strides[3]); + auto pres = reinterpret_cast(data[4] + i * strides[4]); + auto moles = reinterpret_cast(data[5] + i * strides[5]); + auto log_k = reinterpret_cast(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, diff --git a/src/equilibrium/equilibrium_dispatch.hpp b/src/equilibrium/equilibrium_dispatch.hpp index 30423ede..4f7ad1d8 100644 --- a/src/equilibrium/equilibrium_dispatch.hpp +++ b/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); diff --git a/src/equilibrium/equilibrium_options.cpp b/src/equilibrium/equilibrium_options.cpp index a2a7bd33..f3e50efb 100644 --- a/src/equilibrium/equilibrium_options.cpp +++ b/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 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(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(); 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> compositions; std::set 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(); 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(); compositions[name][element] = entry.second.as(); 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 reactions; - for (auto const &reaction_node : config["reactions"]) { + for (auto const& reaction_node : config["reactions"]) { if (reaction_node["type"] && reaction_node["type"].as() != "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; diff --git a/src/equilibrium/phase_equilibrate_tp.h b/src/equilibrium/phase_equilibrate_tp.h index f2ba672c..23b7ac99 100644 --- a/src/equilibrium/phase_equilibrate_tp.h +++ b/src/equilibrium/phase_equilibrate_tp.h @@ -1,22 +1,24 @@ #pragma once #include -#include +#include +#include #include #include #include #include +#include namespace kintera { template -DISPATCH_MACRO T equilibrium_max_error(T const *moles, T pres, - T standard_pressure, T const *log_k, - T const *stoich, int const *phase_ids, +DISPATCH_MACRO T equilibrium_max_error(T const* moles, T pres, + T standard_pressure, T const* log_k, + T const* stoich, int const* phase_ids, int nspecies, int nreaction, int nphase, - int gas_phase, T *phase_totals, - T *residual) { + int gas_phase, T* phase_totals, + T* residual) { for (int p = 0; p < nphase; ++p) phase_totals[p] = 0.; for (int i = 0; i < nspecies; ++i) { phase_totals[phase_ids[i]] += moles[i]; @@ -43,10 +45,10 @@ DISPATCH_MACRO T equilibrium_max_error(T const *moles, T pres, template DISPATCH_MACRO int phase_equilibrate_tp( - T *gain, T *diag, T *out_moles, T temp, T pres, T const *in_moles, - T const *log_k, T const *stoich, int const *phase_ids, int nspecies, + T* gain, T* diag, T* out_moles, T temp, T pres, T const* in_moles, + T const* log_k, T const* stoich, int const* phase_ids, int nspecies, int nreaction, int nphase, int gas_phase, T standard_pressure, T ftol, - T mole_floor, int max_iter, char *work = nullptr) { + T mole_floor, int max_iter, char* work = nullptr) { if (!(temp > 0.) || !(pres > 0.) || nspecies <= 0 || nreaction <= 0 || nphase <= 0 || gas_phase < 0 || gas_phase >= nphase) { diag[0] = 1.; @@ -59,18 +61,18 @@ DISPATCH_MACRO int phase_equilibrate_tp( } } - T *phase_totals, *phase_stoich, *residual, *jac, *constraints, *bounds, *step; - T *trial; + T *phase_totals, *phase_stoich, *residual, *jac, *constraints; + T *bounds, *step, *trial; bool own_work = work == nullptr; if (own_work) { - phase_totals = (T *)malloc(nphase * sizeof(T)); - phase_stoich = (T *)malloc(nphase * nreaction * sizeof(T)); - residual = (T *)malloc(nreaction * sizeof(T)); - jac = (T *)malloc(nreaction * nreaction * sizeof(T)); - constraints = (T *)malloc(nspecies * nreaction * sizeof(T)); - bounds = (T *)malloc(nspecies * sizeof(T)); - step = (T *)malloc(nreaction * sizeof(T)); - trial = (T *)malloc(nspecies * sizeof(T)); + phase_totals = (T*)malloc(nphase * sizeof(T)); + phase_stoich = (T*)malloc(nphase * nreaction * sizeof(T)); + residual = (T*)malloc(nreaction * sizeof(T)); + jac = (T*)malloc(nreaction * nreaction * sizeof(T)); + constraints = (T*)malloc(nspecies * nreaction * sizeof(T)); + bounds = (T*)malloc(nspecies * sizeof(T)); + step = (T*)malloc(nreaction * sizeof(T)); + trial = (T*)malloc(nspecies * sizeof(T)); } else { phase_totals = alloc_from(work, nphase); phase_stoich = alloc_from(work, nphase * nreaction); @@ -90,7 +92,7 @@ DISPATCH_MACRO int phase_equilibrate_tp( phase_stoich[phase_ids[i] * nreaction + j] += stoich[i * nreaction + j]; } } - T target = log(1. + ftol); + T target = fmax(log1p(ftol), 4. * std::numeric_limits::epsilon()); T max_error = 0.; int status = 2; int iter = 0; @@ -127,26 +129,22 @@ DISPATCH_MACRO int phase_equilibrate_tp( } } + // Use the original square Jacobian whenever its direction can be made + // feasible by line search. KKT remains the active-bound fallback. int kkt_iter = max_iter; - int err = leastsq_kkt(step, jac, constraints, bounds, nreaction, nreaction, - nspecies, 0, &kkt_iter, 1.e-12, work); + int err = constrained_newton_step(step, jac, constraints, bounds, nreaction, + nspecies, &kkt_iter, 1.e-12, work); if (err != 0) { status = 3; break; } - T scale = 1.; bool accepted = false; - while (scale >= 1.e-8) { - bool positive = true; - for (int i = 0; i < nspecies; ++i) { - T delta = 0.; - for (int j = 0; j < nreaction; ++j) { - delta += stoich[i * nreaction + j] * step[j]; - } - trial[i] = out_moles[i] + scale * delta; - if (!(trial[i] > mole_floor)) positive = false; - } + T scale = 1.; + while (scale >= 1.e-12) { + bool positive = constrained_newton_trial(trial, out_moles, constraints, + step, nspecies, nreaction, + nspecies, 0, scale, mole_floor); if (positive) { T trial_error = equilibrium_max_error( trial, pres, standard_pressure, log_k, stoich, phase_ids, nspecies, diff --git a/src/math/constrained_newton.h b/src/math/constrained_newton.h new file mode 100644 index 00000000..cbb1c1cf --- /dev/null +++ b/src/math/constrained_newton.h @@ -0,0 +1,98 @@ +#pragma once + +// C/C++ +#include +#include +#include + +// base +#include + +// kintera +#include +#include +#include + +namespace kintera { + +/*! + * \brief Solve a square Newton system, retaining KKT for active bounds. + * + * The direct solve avoids the loss of precision caused by forming A^T.A. A + * direction that violates a bound at unit length is still usable when a + * positive scaling makes it feasible. Only a direction pointing out of an + * already-active bound is sent through constrained least squares. + */ +template +DISPATCH_MACRO int constrained_newton_step(T* b, T const* a, T const* c, + T const* d, int n, int nconstraint, + int* max_iter, float reg = 0., + char* work = nullptr) { + T *direct_a, *direct_b; + if (work == nullptr) { + direct_a = (T*)malloc(n * n * sizeof(T)); + direct_b = (T*)malloc(n * sizeof(T)); + } else { + char* cursor = work; + direct_a = alloc_from(cursor, n * n); + direct_b = alloc_from(cursor, n); + } + + memcpy(direct_a, a, n * n * sizeof(T)); + memcpy(direct_b, b, n * sizeof(T)); + bool usable = dsolve_lu(direct_a, direct_b, n); + for (int j = 0; j < n && usable; ++j) { + if (!std::isfinite(direct_b[j])) usable = false; + } + for (int i = 0; i < nconstraint && usable; ++i) { + T value = 0.; + for (int j = 0; j < n; ++j) value += c[i * n + j] * direct_b[j]; + if (!(d[i] > 0.) && value > 0.) usable = false; + } + + if (usable) { + memcpy(b, direct_b, n * sizeof(T)); + *max_iter = 1; + } + if (work == nullptr) { + free(direct_a); + free(direct_b); + } + if (usable) return 0; + + return leastsq_kkt(b, a, c, d, n, n, nconstraint, 0, max_iter, reg, work); +} + +/*! + * \brief Form a scaled Newton trial state and enforce simple bounds. + * + * `constraint` uses the same convention as leastsq_kkt, so the state update + * is state - constraint.step. The first `npositive` entries must remain + * strictly above `floor`; all remaining entries must remain non-negative. + * A positive `max_ratio` additionally limits relative changes of positive + * entries in the first `nlimited` entries. + */ +template +DISPATCH_MACRO bool constrained_newton_trial(T* trial, T const* state, + T const* constraint, T const* step, + int nstate, int nstep, + int npositive, int nlimited, + T scale, T floor = 0., + T max_ratio = 0.) { + for (int i = 0; i < nstate; ++i) { + T delta = 0.; + for (int j = 0; j < nstep; ++j) { + delta -= constraint[i * nstep + j] * step[j]; + } + trial[i] = state[i] + scale * delta; + if (!std::isfinite(trial[i]) || trial[i] < 0.) return false; + if (i < npositive && !(trial[i] > floor)) return false; + if (max_ratio > 0. && i < nlimited && state[i] > 0. && + (trial[i] > max_ratio * state[i] || trial[i] * max_ratio < state[i])) { + return false; + } + } + return true; +} + +} // namespace kintera diff --git a/src/species.cpp b/src/species.cpp index 62415768..45a66aa0 100644 --- a/src/species.cpp +++ b/src/species.cpp @@ -55,14 +55,14 @@ void clear_species_registry() { } // namespace -static std::unordered_map &get_nasa9_db() { +static std::unordered_map& get_nasa9_db() { static std::unordered_map db; static std::once_flag initialized; std::call_once(initialized, [&] { std::string path; try { path = find_resource("nasa9.dat"); - } catch (std::exception const &e) { + } catch (std::exception const& e) { TORCH_CHECK(false, e.what()); } std::ifstream ifs(path); @@ -100,15 +100,15 @@ static std::unordered_map &get_nasa9_db() { } at::Tensor nasa9_gibbs_rt(at::Tensor temp, - std::vector const &species) { + std::vector const& species) { TORCH_CHECK(temp.is_floating_point(), "temp must be a floating-point tensor"); TORCH_CHECK(!species.empty(), "NASA-9 species list must not be empty"); - auto const &database = get_nasa9_db(); + auto const& database = get_nasa9_db(); Nasa9CoeffTable low; Nasa9CoeffTable high; low.reserve(species.size()); high.reserve(species.size()); - for (auto const &name : species) { + for (auto const& name : species) { auto found = database.find(name); TORCH_CHECK(found != database.end(), "NASA-9 species not found: ", name); low.push_back(found->second.low); @@ -134,18 +134,18 @@ void init_species_from_yaml(std::string filename) { init_species_from_yaml(config); } -void init_species_from_yaml(YAML::Node const &config) { +void init_species_from_yaml(YAML::Node const& config) { // check if species are defined TORCH_CHECK(config["species"], "'species' is not defined in the kintera configuration file"); clear_species_registry(); - for (const auto &sp : config["species"]) { + for (const auto& sp : config["species"]) { species_names.push_back(sp["name"].as()); std::map comp; - for (const auto &it : sp["composition"]) { + for (const auto& it : sp["composition"]) { std::string key = it.first.as(); double value = it.second.as(); comp[key] = value; @@ -175,7 +175,7 @@ void init_species_from_yaml(YAML::Node const &config) { std::array high_coeffs = {}; double Tmid = 1000.0; - auto &nasa9_db = get_nasa9_db(); + auto& nasa9_db = get_nasa9_db(); auto name = sp["name"].as(); auto it = nasa9_db.find(name); if (it != nasa9_db.end()) { @@ -193,13 +193,13 @@ void init_species_from_yaml(YAML::Node const &config) { species_initialized = true; } -void ensure_species_initialized(std::string const &filename) { +void ensure_species_initialized(std::string const& filename) { if (!species_initialized) { init_species_from_yaml(filename); } } -void ensure_species_initialized(YAML::Node const &config) { +void ensure_species_initialized(YAML::Node const& config) { if (!species_initialized) { init_species_from_yaml(config); } @@ -222,7 +222,7 @@ std::vector SpeciesThermoImpl::species() const { } at::Tensor SpeciesThermoImpl::narrow_copy(at::Tensor data, - SpeciesThermo const &other) const { + SpeciesThermo const& other) const { auto source_ids = merge_vectors(vapor_ids(), cloud_ids()); auto other_ids = merge_vectors(other->vapor_ids(), other->cloud_ids()); std::vector indices; @@ -241,9 +241,9 @@ at::Tensor SpeciesThermoImpl::narrow_copy(at::Tensor data, return data.index_select(-1, id); } -void SpeciesThermoImpl::accumulate(at::Tensor &data, - at::Tensor const &other_data, - SpeciesThermo const &other) const { +void SpeciesThermoImpl::accumulate(at::Tensor& data, + at::Tensor const& other_data, + SpeciesThermo const& other) const { auto source_ids = merge_vectors(vapor_ids(), cloud_ids()); auto other_ids = merge_vectors(other->vapor_ids(), other->cloud_ids()); std::vector indices; @@ -262,12 +262,12 @@ void SpeciesThermoImpl::accumulate(at::Tensor &data, } bool SpeciesThermoImpl::has_nasa9() const { - for (auto const &coeffs : nasa9_low()) { + for (auto const& coeffs : nasa9_low()) { for (double v : coeffs) { if (v != 0.0) return true; } } - for (auto const &coeffs : nasa9_high()) { + for (auto const& coeffs : nasa9_high()) { for (double v : coeffs) { if (v != 0.0) return true; } @@ -276,8 +276,8 @@ bool SpeciesThermoImpl::has_nasa9() const { } static at::Tensor nasa9_coeffs_to_tensor( - std::vector> const &coeffs, - c10::TensorOptions const &options) { + std::vector> const& coeffs, + c10::TensorOptions const& options) { auto tensor = torch::empty({static_cast(coeffs.size()), 9}, torch::dtype(torch::kFloat64)); if (!coeffs.empty()) { @@ -292,17 +292,17 @@ static at::Tensor nasa9_coeffs_to_tensor( } at::Tensor SpeciesThermoImpl::nasa9_coeffs_low_tensor( - c10::TensorOptions const &options) const { + c10::TensorOptions const& options) const { return nasa9_coeffs_to_tensor(nasa9_low(), options); } at::Tensor SpeciesThermoImpl::nasa9_coeffs_high_tensor( - c10::TensorOptions const &options) const { + c10::TensorOptions const& options) const { return nasa9_coeffs_to_tensor(nasa9_high(), options); } at::Tensor SpeciesThermoImpl::nasa9_Tmid_tensor( - c10::TensorOptions const &options) const { + c10::TensorOptions const& options) const { auto tensor = torch::empty({static_cast(nasa9_Tmid().size())}, torch::dtype(torch::kFloat64)); if (!nasa9_Tmid().empty()) { @@ -353,7 +353,7 @@ void populate_thermo(SpeciesThermo thermo) { } } -void check_dimensions(SpeciesThermo const &thermo) { +void check_dimensions(SpeciesThermo const& thermo) { int nspecies = thermo->vapor_ids().size() + thermo->cloud_ids().size(); TORCH_CHECK(thermo->cref_R().size() == nspecies, @@ -404,8 +404,8 @@ void check_dimensions(SpeciesThermo const &thermo) { ". Expected = ", nspecies); } -SpeciesThermo merge_thermo(SpeciesThermo const &thermo1, - SpeciesThermo const &thermo2) { +SpeciesThermo merge_thermo(SpeciesThermo const& thermo1, + SpeciesThermo const& thermo2) { // check dimensions check_dimensions(thermo1); check_dimensions(thermo2); @@ -413,20 +413,20 @@ SpeciesThermo merge_thermo(SpeciesThermo const &thermo1, // return a new SpeciesThermo object with merged data auto merged = SpeciesThermoImpl::create(); - auto &vapor_ids = merged->vapor_ids(); - auto &cloud_ids = merged->cloud_ids(); - - auto &cref_R = merged->cref_R(); - auto &uref_R = merged->uref_R(); - auto &sref_R = merged->sref_R(); - auto &intEng_R_extra = merged->intEng_R_extra(); - auto &cp_R_extra = merged->cp_R_extra(); - auto &entropy_R_extra = merged->entropy_R_extra(); - auto &czh = merged->czh(); - auto &czh_ddC = merged->czh_ddC(); - auto &nasa9_low = merged->nasa9_low(); - auto &nasa9_high = merged->nasa9_high(); - auto &nasa9_Tmid = merged->nasa9_Tmid(); + auto& vapor_ids = merged->vapor_ids(); + auto& cloud_ids = merged->cloud_ids(); + + auto& cref_R = merged->cref_R(); + auto& uref_R = merged->uref_R(); + auto& sref_R = merged->sref_R(); + auto& intEng_R_extra = merged->intEng_R_extra(); + auto& cp_R_extra = merged->cp_R_extra(); + auto& entropy_R_extra = merged->entropy_R_extra(); + auto& czh = merged->czh(); + auto& czh_ddC = merged->czh_ddC(); + auto& nasa9_low = merged->nasa9_low(); + auto& nasa9_high = merged->nasa9_high(); + auto& nasa9_Tmid = merged->nasa9_Tmid(); // concatenate fields int nvapor1 = thermo1->vapor_ids().size(); @@ -534,7 +534,7 @@ SpeciesThermo merge_thermo(SpeciesThermo const &thermo1, cloud_ids = sort_vectors(cloud_ids, cidx); // add nvapor to cidx - for (auto &idx : cidx) idx += nvapor; + for (auto& idx : cidx) idx += nvapor; auto sorted = merge_vectors(vidx, cidx); diff --git a/src/species.hpp b/src/species.hpp index 00bf8b64..96234d2e 100644 --- a/src/species.hpp +++ b/src/species.hpp @@ -28,15 +28,15 @@ using Nasa9CoeffArray = std::array; using Nasa9CoeffTable = std::vector; void init_species_from_yaml(std::string filename); -void init_species_from_yaml(YAML::Node const &config); +void init_species_from_yaml(YAML::Node const& config); //! Initialize species and thermo data from a KINETICS-base master input file. -void init_species_from_kinetics_base(std::string const &master_input_path); -void ensure_species_initialized(std::string const &filename); -void ensure_species_initialized(YAML::Node const &config); +void init_species_from_kinetics_base(std::string const& master_input_path); +void ensure_species_initialized(std::string const& filename); +void ensure_species_initialized(YAML::Node const& config); //! Evaluate standard-state Gibbs energy g/RT from the bundled NASA-9 database. at::Tensor nasa9_gibbs_rt(at::Tensor temp, - std::vector const &species); + std::vector const& species); struct SpeciesThermoImpl { static std::shared_ptr create() { @@ -49,16 +49,16 @@ struct SpeciesThermoImpl { std::vector species() const; at::Tensor narrow_copy(at::Tensor data, - std::shared_ptr const &other) const; - void accumulate(at::Tensor &data, at::Tensor const &other_data, - std::shared_ptr const &other) const; + std::shared_ptr const& other) const; + void accumulate(at::Tensor& data, at::Tensor const& other_data, + std::shared_ptr const& other) const; bool has_nasa9() const; at::Tensor nasa9_coeffs_low_tensor( - c10::TensorOptions const &options = c10::TensorOptions()) const; + c10::TensorOptions const& options = c10::TensorOptions()) const; at::Tensor nasa9_coeffs_high_tensor( - c10::TensorOptions const &options = c10::TensorOptions()) const; + c10::TensorOptions const& options = c10::TensorOptions()) const; at::Tensor nasa9_Tmid_tensor( - c10::TensorOptions const &options = c10::TensorOptions()) const; + c10::TensorOptions const& options = c10::TensorOptions()) const; ADD_ARG(std::vector, vapor_ids); ADD_ARG(std::vector, cloud_ids); @@ -96,10 +96,10 @@ using SpeciesThermo = std::shared_ptr; void populate_thermo(SpeciesThermo thermo); -void check_dimensions(SpeciesThermo const &thermo); +void check_dimensions(SpeciesThermo const& thermo); -SpeciesThermo merge_thermo(SpeciesThermo const &thermo1, - SpeciesThermo const &thermo2); +SpeciesThermo merge_thermo(SpeciesThermo const& thermo1, + SpeciesThermo const& thermo2); extern std::vector species_names; extern std::vector species_weights; diff --git a/src/thermo/equilibrate_tp.h b/src/thermo/equilibrate_tp.h index 28b2cfd6..1b1cf08f 100644 --- a/src/thermo/equilibrate_tp.h +++ b/src/thermo/equilibrate_tp.h @@ -10,8 +10,8 @@ #include // kintera +#include #include -#include #include @@ -145,7 +145,6 @@ DISPATCH_MACRO int equilibrate_tp(T* gain, T* diag, T* xfrac, T temp, T pres, int iter = 0; int kkt_err = 0; - T lambda = 0.; // rate scale factor while (iter++ < *max_iter) { /*printf("iter = %d\n ", iter); // print xfrac @@ -242,10 +241,16 @@ DISPATCH_MACRO int equilibrate_tp(T* gain, T* diag, T* xfrac, T temp, T pres, } // note that stoich_active is negated - // solve constrained optimization problem (KKT) + T current_error = 0.; + for (int k = 0; k < (*nactive); ++k) { + current_error = fmax(current_error, fabs(rhs[k])); + } + + // Solve the square Newton system directly unless an active bound requires + // the constrained KKT fallback. int max_kkt_iter = *max_iter; - kkt_err = leastsq_kkt(rhs, gain, stoich_active, xfrac, *nactive, *nactive, - nspecies, 0, &max_kkt_iter, 0., work); + kkt_err = constrained_newton_step(rhs, gain, stoich_active, xfrac, *nactive, + nspecies, &max_kkt_iter, 0., work); if (kkt_err != 0) break; /* print rate @@ -255,28 +260,50 @@ DISPATCH_MACRO int equilibrate_tp(T* gain, T* diag, T* xfrac, T temp, T pres, } printf("\n");*/ - // rate -> xfrac - // copy xfrac to xfrac0 - memcpy(xfrac0, xfrac, nspecies * sizeof(T)); + // Backtrack until the bounded trial state reduces the nonlinear + // equilibrium residual. Normalization does not change gas mixing ratios. + bool accepted = false; T lambda = 1.; // scale - T xsum; - while (true) { - bool positive_vapor = true; - xsum = 0.; - for (int i = 0; i < nspecies; i++) { - for (int k = 0; k < (*nactive); k++) { - xfrac[i] -= stoich_active[i * (*nactive) + k] * rhs[k] * lambda; + while (lambda >= 1.e-12) { + bool feasible = + constrained_newton_trial(xfrac0, xfrac, stoich_active, rhs, nspecies, + *nactive, ngas, 0, lambda); + T xsum = 0.; + for (int i = 0; i < nspecies && feasible; ++i) xsum += xfrac0[i]; + feasible = feasible && xsum > 0. && std::isfinite(xsum); + for (int i = 0; i < nspecies && feasible; ++i) xfrac0[i] /= xsum; + + T trial_error = 0.; + T trial_xg = 0.; + for (int i = 0; i < ngas && feasible; ++i) trial_xg += xfrac0[i]; + feasible = feasible && trial_xg > 0.; + for (int k = 0; k < (*nactive) && feasible; ++k) { + int j = reaction_set[k]; + T log_frac_sum = 0.; + for (int i = 0; i < nspecies; ++i) { + if (stoich[i * nreaction + j] < 0.) { + if (!(xfrac0[i] > 0.)) { + feasible = false; + break; + } + log_frac_sum += + (-stoich[i * nreaction + j]) * log(xfrac0[i] / trial_xg); + } } - if (i < ngas && xfrac[i] <= 0.) positive_vapor = false; - xsum += xfrac[i]; + if (!feasible) break; + trial_error = fmax(trial_error, fabs(logsvp[j] - log_frac_sum)); } - if (positive_vapor) break; - lambda *= 0.99; - memcpy(xfrac, xfrac0, nspecies * sizeof(T)); + if (feasible && trial_error < current_error) { + memcpy(xfrac, xfrac0, nspecies * sizeof(T)); + accepted = true; + break; + } + lambda *= .5; + } + if (!accepted) { + kkt_err = 3; + break; } - - // re-normalize mole fractions - for (int i = 0; i < nspecies; i++) xfrac[i] /= xsum; } /*///////// Construct a gain matrix of active reactions /////////// diff --git a/src/thermo/equilibrate_uv.h b/src/thermo/equilibrate_uv.h index 7610b1fb..bfa0a105 100644 --- a/src/thermo/equilibrate_uv.h +++ b/src/thermo/equilibrate_uv.h @@ -10,8 +10,8 @@ // kintera #include +#include #include -#include #include @@ -240,59 +240,101 @@ DISPATCH_MACRO int equilibrate_uv( } // note that stoich_active is negated - // solve constrained optimization problem (KKT) + T current_error = 0.; + for (int k = 0; k < (*nactive); ++k) { + current_error = fmax(current_error, fabs(rhs[k])); + } + + // Solve the original square Newton system, using KKT only when a bound is + // active in the Newton direction. int max_kkt_iter = *max_iter; - err_code = leastsq_kkt(rhs, gain, stoich_active, conc, *nactive, *nactive, - nspecies, 0, &max_kkt_iter, -1.e-10, work); + err_code = constrained_newton_step(rhs, gain, stoich_active, conc, *nactive, + nspecies, &max_kkt_iter, -1.e-10, work); if (err_code != 0) break; - // rate -> conc - memcpy(conc0, conc, nspecies * sizeof(T)); + // Backtrack on the coupled UV residual. Each trial composition gets its + // own energy-conserving temperature before its chemistry is evaluated. + bool accepted = false; T lambda = 1.; // scale - while (true) { - bool good = true; - for (int i = 0; i < nspecies; i++) { - for (int k = 0; k < (*nactive); k++) { - conc[i] -= stoich_active[i * (*nactive) + k] * rhs[k] * lambda; + T accepted_temp = *temp; + while (lambda >= 1.e-12) { + bool feasible = constrained_newton_trial( + conc0, conc, stoich_active, rhs, nspecies, *nactive, 0, ngas, lambda, + static_cast(0.), static_cast(100.)); + T trial_temp = *temp; + bool energy_converged = false; + for (int temp_iter = 0; temp_iter < 50 && feasible; ++temp_iter) { + T zh = 0.; + T zc = 0.; + for (int i = 0; i < nspecies; ++i) { + T energy = intEng_offset[i] + cv_const[i] * trial_temp; + if (intEng_R_extra[i]) { + energy += intEng_R_extra[i](trial_temp, conc0[i]) * constants::Rgas; + } + T heat = cv_const[i]; + if (cv_R_extra[i]) { + heat += cv_R_extra[i](trial_temp, conc0[i]) * constants::Rgas; + } + zh += energy * conc0[i]; + zc += heat * conc0[i]; + } + feasible = zc > 0. && std::isfinite(zc); + if (!feasible) break; + T delta_temp = (h0 - zh) / zc; + trial_temp += delta_temp; + feasible = trial_temp > 0. && std::isfinite(trial_temp); + if (fabs(delta_temp) <= 1.e-4) { + energy_converged = feasible; + break; } - if ((i < ngas) && (conc0[i] > 0.) && - ((conc[i] / conc0[i] > 100.) || (conc[i] / conc0[i] < 0.01))) - good = false; } - if (good) break; - lambda *= 0.99; - memcpy(conc, conc0, nspecies * sizeof(T)); - } - - // temperature iteration - T temp0 = 0.; - while (fabs(*temp - temp0) > 1e-4) { - T zh = 0.; - T zc = 0.; + feasible = feasible && energy_converged; - // re-evaluate internal energy and its derivative - for (int i = 0; i < nspecies; i++) { - intEng[i] = intEng_offset[i] + cv_const[i] * (*temp); - if (intEng_R_extra[i]) { - intEng[i] += intEng_R_extra[i](*temp, conc[i]) * constants::Rgas; - } - intEng_ddT[i] = cv_const[i]; - if (cv_R_extra[i]) { - intEng_ddT[i] += cv_R_extra[i](*temp, conc[i]) * constants::Rgas; + T trial_error = 0.; + for (int k = 0; k < (*nactive) && feasible; ++k) { + int j = reaction_set[k]; + T stoich_sum = 0.; + T log_conc_sum = 0.; + for (int i = 0; i < nspecies; ++i) { + if (stoich[i * nreaction + j] < 0.) { + stoich_sum += -stoich[i * nreaction + j]; + if (conc0[i] == 0.) { + log_conc_sum = -99.; + } else { + log_conc_sum += (-stoich[i * nreaction + j]) * log(conc0[i]); + } + } } - zh += intEng[i] * conc[i]; - zc += intEng_ddT[i] * conc[i]; + T trial_logsvp = logsvp_func[j](trial_temp) - + stoich_sum * log(constants::Rgas * trial_temp); + trial_error = fmax(trial_error, fabs(trial_logsvp - log_conc_sum)); } - - temp0 = *temp; - (*temp) += (h0 - zh) / zc; + if (feasible && trial_error < current_error) { + memcpy(conc, conc0, nspecies * sizeof(T)); + accepted_temp = trial_temp; + accepted = true; + break; + } + lambda *= .5; } - - if (*temp <= 0.) { - printf("Error: Non-positive temperature after adjustment.\n"); - err_code = 3; // error: non-positive temperature after adjustment + if (!accepted) { + err_code = 4; break; } + *temp = accepted_temp; + + // Keep the thermodynamic linearization synchronized with the accepted + // temperature and composition for the next coupled Newton iteration. + for (int i = 0; i < nspecies; ++i) { + intEng[i] = intEng_offset[i] + cv_const[i] * (*temp); + if (intEng_R_extra[i]) { + intEng[i] += intEng_R_extra[i](*temp, conc[i]) * constants::Rgas; + } + intEng_ddT[i] = cv_const[i]; + if (cv_R_extra[i]) { + intEng_ddT[i] += cv_R_extra[i](*temp, conc[i]) * constants::Rgas; + } + } } // restore the reaction order of gain diff --git a/src/utils/molar_mass.cpp b/src/utils/molar_mass.cpp index d096d3b4..bd8a28b2 100644 --- a/src/utils/molar_mass.cpp +++ b/src/utils/molar_mass.cpp @@ -10,23 +10,23 @@ namespace kintera { -double atomic_mass(std::string const &element) { +double atomic_mass(std::string const& element) { return harp::get_element_weight(element) * 1.e-3; } -double molar_mass(Composition const &composition) { +double molar_mass(Composition const& composition) { return harp::get_compound_weight(composition); } std::vector molar_masses( - std::vector const &elements, - std::vector> const &element_matrix) { + std::vector const& elements, + std::vector> const& element_matrix) { TORCH_CHECK(elements.size() == element_matrix.size(), "elements and element_matrix row counts must match"); if (elements.empty()) return {}; auto ncomponent = element_matrix[0].size(); - for (auto const &row : element_matrix) { + for (auto const& row : element_matrix) { TORCH_CHECK(row.size() == ncomponent, "element_matrix rows must have equal length"); } @@ -45,29 +45,29 @@ std::vector molar_masses( return result; } -std::vector molar_masses_from_yaml(std::string const &filename) { +std::vector molar_masses_from_yaml(std::string const& filename) { auto config = YAML::LoadFile(find_resource(filename)); TORCH_CHECK(config["phases"], "chemistry YAML requires 'phases'"); TORCH_CHECK(config["species"], "chemistry YAML requires 'species'"); std::map compositions; - 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(); TORCH_CHECK(compositions.emplace(name, Composition{}).second, "duplicate species definition: ", name); - for (auto const &entry : species_node["composition"]) { + for (auto const& entry : species_node["composition"]) { compositions.at(name)[entry.first.as()] = entry.second.as(); } } std::vector result; - for (auto const &phase_node : config["phases"]) { + for (auto const& phase_node : config["phases"]) { TORCH_CHECK(phase_node["species"], "phase is missing 'species'"); - for (auto const &component_node : phase_node["species"]) { + for (auto const& component_node : phase_node["species"]) { auto component = component_node.as(); auto found = compositions.find(component); TORCH_CHECK(found != compositions.end(), diff --git a/src/utils/molar_mass.hpp b/src/utils/molar_mass.hpp index 6531496c..37c3b6be 100644 --- a/src/utils/molar_mass.hpp +++ b/src/utils/molar_mass.hpp @@ -7,17 +7,17 @@ namespace kintera { //! Atomic molar mass [kg/mol]. -double atomic_mass(std::string const &element); +double atomic_mass(std::string const& element); //! Compound molar mass [kg/mol]. -double molar_mass(Composition const &composition); +double molar_mass(Composition const& composition); //! Component molar masses [kg/mol] from an element-by-component matrix. std::vector molar_masses( - std::vector const &elements, - std::vector> const &element_matrix); + std::vector const& elements, + std::vector> const& element_matrix); //! Component molar masses [kg/mol] in phase order from a chemistry YAML file. -std::vector molar_masses_from_yaml(std::string const &filename); +std::vector molar_masses_from_yaml(std::string const& filename); } // namespace kintera diff --git a/tests/test_equilibrium.cpp b/tests/test_equilibrium.cpp index d32e6618..aeb75a97 100644 --- a/tests/test_equilibrium.cpp +++ b/tests/test_equilibrium.cpp @@ -135,7 +135,7 @@ INSTANTIATE_TEST_SUITE_P( Parameters{torch::kCPU, torch::kFloat64}, Parameters{torch::kCUDA, torch::kFloat32}, Parameters{torch::kCUDA, torch::kFloat64}), - [](const testing::TestParamInfo &info) { + [](const testing::TestParamInfo& info) { std::string name = torch::Device(info.param.device_type).str(); name += "_"; name += torch::toString(info.param.dtype); diff --git a/tests/test_equilibrium.py b/tests/test_equilibrium.py index 1ed6314e..923124e5 100644 --- a/tests/test_equilibrium.py +++ b/tests/test_equilibrium.py @@ -65,3 +65,86 @@ def test_nasa9_yaml_equilibrium_case(): result @ torch.tensor([0.0, 2.0, 1.0], dtype=result.dtype), initial @ torch.tensor([0.0, 2.0, 1.0], dtype=initial.dtype), ) + + +def test_ill_conditioned_multiphase_newton_step(): + """Regression for normal-equation precision loss in a Young22 state.""" + components = [ + "MgO_silicate", "SiO2_silicate", "MgSiO3_silicate", + "FeO_silicate", "FeSiO3_silicate", "Na2O_silicate", + "Na2SiO3_silicate", "H2_silicate", "H2O_silicate", + "CO_silicate", "CO2_silicate", "Fe_metal", "Si_metal", + "O_metal", "H_metal", "H2_gas", "CO_gas", "CO2_gas", + "CH4_gas", "O2_gas", "H2O_gas", "Fe_gas", "Mg_gas", + "SiO_gas", "Na_gas", + ] + reactions = [ + "Na2SiO3_silicate <=> Na2O_silicate + SiO2_silicate", + "0.5 SiO2_silicate + Fe_metal <=> FeO_silicate + 0.5 Si_metal", + "MgSiO3_silicate <=> MgO_silicate + SiO2_silicate", + "O_metal + 0.5 Si_metal <=> 0.5 SiO2_silicate", + "2 H_metal <=> H2_silicate", + "FeSiO3_silicate <=> FeO_silicate + SiO2_silicate", + "2 H2O_silicate + Si_metal <=> SiO2_silicate + 2 H2_silicate", + "CO_gas + 0.5 O2_gas <=> CO2_gas", + "CH4_gas + 0.5 O2_gas <=> 2 H2_gas + CO_gas", + "H2_gas + 0.5 O2_gas <=> H2O_gas", + "FeO_silicate <=> Fe_gas + 0.5 O2_gas", + "MgO_silicate <=> Mg_gas + 0.5 O2_gas", + "SiO2_silicate <=> SiO_gas + 0.5 O2_gas", + "Na2O_silicate <=> 2 Na_gas + 0.5 O2_gas", + "H2_gas <=> H2_silicate", + "H2O_gas <=> H2O_silicate", + "CO_gas <=> CO_silicate", + "CO2_gas <=> CO2_silicate", + ] + options = ( + EquilibriumOptions() + .components(components) + .phases(["silicate", "metal", "gas"]) + .phase_ids([0] * 11 + [1] * 4 + [2] * 10) + .reactions(reactions) + .gas_phase(2) + .standard_pressure(1.0e5) + .max_iter(1000) + .ftol(1.0e-8) + .mole_floor(1.0e-40) + ) + initial = torch.tensor( + [ + 229.45501199233973, 7.170469124760617, 6604.002063904528, + 3.5852345623803084, 250.96641936662164, 50.19328387332432, + 7.170469124760617, 7.1704691247606165e-9, + 7.1704691247606165e-9, 7.1704691247606165e-9, + 7.1704691247606165e-9, 4476.516196031982, + 4.476516196031981e-9, 4.476516196031981e-9, + 4.476516196031981e-9, 19570.252531955153, 19.589842374329482, + 1.9589842374329483e-8, 1.9589842374329483e-8, + 1.9589842374329483e-8, 1.9589842374329483e-8, + 1.9589842374329483e-8, 1.9589842374329483e-8, + 1.9589842374329483e-8, 1.9589842374329483e-8, + ], + dtype=torch.float64, + ) + log_k = torch.tensor( + [ + -4.034640768501788, -7.62490310519855, -2.039313234293211, + 7.218657926794565, 0.7959817637548838, -0.1371317344272009, + 26.01075006664782, -4.50470154179624, 38.86385560946831, + -3.708764843086439, 4.050968036542196, 1.772570237146461, + 5.579140721827175, 14.1841721338006, -13.701917787416148, + -14.27932134370162, -16.21258232342021, -14.87214162753718, + ], + dtype=torch.float64, + ) + result, _, diagnostics = EquilibriumTP(options)( + torch.tensor(4500.0, dtype=torch.float64), + torch.tensor(12599210499.047432, dtype=torch.float64), + initial, + log_k, + ) + + assert diagnostics[0] == 0 + assert diagnostics[2] < 1.0e-8 + assert torch.isfinite(result).all() + assert torch.all(result > 0.0) diff --git a/tests/test_thermo.cpp b/tests/test_thermo.cpp index c1d5f45c..d8795ac0 100644 --- a/tests/test_thermo.cpp +++ b/tests/test_thermo.cpp @@ -6,6 +6,7 @@ // kintera #include +#include #include #include @@ -20,6 +21,44 @@ using namespace kintera; using namespace torch::indexing; +TEST(LeastSquaresKkt, DirectlySolvesIllConditionedSquareSystem) { + double matrix[] = {1., 1., 1., 1. + 1.e-10}; + double rhs[] = {0., -1.e-10}; + int max_iter = 10; + + int status = constrained_newton_step( + rhs, matrix, static_cast(nullptr), + static_cast(nullptr), 2, 0, &max_iter); + + EXPECT_EQ(status, 0); + EXPECT_NEAR(rhs[0], 1., 1.e-6); + EXPECT_NEAR(rhs[1], -1., 1.e-6); + EXPECT_EQ(max_iter, 1); +} + +TEST(LeastSquaresKkt, LineSearchCanScaleNewtonDirectionToBound) { + double matrix[] = {1.}; + double constraint[] = {1.}; + double bound[] = {1.}; + double rhs[] = {2.}; + int max_iter = 10; + + int status = + constrained_newton_step(rhs, matrix, constraint, bound, 1, 1, &max_iter); + + EXPECT_EQ(status, 0); + EXPECT_DOUBLE_EQ(rhs[0], 2.); + EXPECT_EQ(max_iter, 1); + + double state[] = {1.}; + double trial[] = {0.}; + EXPECT_FALSE( + constrained_newton_trial(trial, state, constraint, rhs, 1, 1, 1, 0, 1.)); + EXPECT_TRUE(constrained_newton_trial(trial, state, constraint, rhs, 1, 1, 1, + 0, 0.25)); + EXPECT_DOUBLE_EQ(trial[0], 0.5); +} + TEST_P(DeviceTest, thermo_y) { auto op_thermo = ThermoOptionsImpl::from_yaml("jupiter.yaml");