diff --git a/poet/include/EquilibriumWrapper.hpp b/poet/include/EquilibriumWrapper.hpp new file mode 100644 index 00000000..a22c1cd1 --- /dev/null +++ b/poet/include/EquilibriumWrapper.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include "PPassemblage.h" +#include "WrapperBase.hpp" +#include + +class EquilibriumWrapper : public WrapperBase { +public: + EquilibriumWrapper(cxxPPassemblage *ppassemblage, + const std::vector &ppassemblage_names); + + void get(std::span &data) const override; + + void set(const std::span &data) override; + + static std::vector + names(const cxxPPassemblage *ppassemblage, + std::vector &ppassemblage_names); + +private: + cxxPPassemblage *ppassemblage; + + class EquilibriumCompWrapper : public WrapperBase { + public: + EquilibriumCompWrapper(cxxPPassemblageComp &comp); + + void get(std::span &data) const override; + + void set(const std::span &data) override; + + static std::vector names(const cxxPPassemblageComp &comp); + + private: + cxxPPassemblageComp ∁ + }; + + std::vector> equilibrium_comps; +}; \ No newline at end of file diff --git a/poet/include/ExchangeWrapper.hpp b/poet/include/ExchangeWrapper.hpp new file mode 100644 index 00000000..dd4824cf --- /dev/null +++ b/poet/include/ExchangeWrapper.hpp @@ -0,0 +1,43 @@ +#pragma once + +#include "ExchComp.h" +#include "Exchange.h" +#include "WrapperBase.hpp" +#include "phrqtype.h" +#include +#include +#include +#include + +#include + +class ExchangeWrapper : public WrapperBase { +public: + ExchangeWrapper(cxxExchange *exch, const std::vector &exchanger); + + void get(std::span &exchange) const; + + void set(const std::span &exchange); + + static std::vector + names(cxxExchange *exchange, std::vector &exchange_formulas); + +private: + cxxExchange *exchange; + + class ExchangeCompWrapper : public WrapperBase { + public: + ExchangeCompWrapper(cxxExchComp &comp); + + void get(std::span &exchange) const; + + void set(const std::span &exchange); + + private: + cxxExchComp &exch_comp; + + static constexpr std::size_t NUM_NOT_TOTALS = 5; + }; + + std::vector> exchange_comps; +}; \ No newline at end of file diff --git a/poet/include/IPhreeqcPOET.hpp b/poet/include/IPhreeqcPOET.hpp deleted file mode 100644 index bb3812bb..00000000 --- a/poet/include/IPhreeqcPOET.hpp +++ /dev/null @@ -1,192 +0,0 @@ -#pragma once - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -enum { POET_SOL = 0, POET_EXCH, POET_KIN, POET_EQUIL, POET_SURF }; - -class IPhreeqcPOET : public IPhreeqc { -public: - IPhreeqcPOET(const std::string &database, const std::string &input_script) - : IPhreeqc() { - this->LoadDatabaseString(database.c_str()); - this->RunString(input_script.c_str()); - - this->parseInit(); - } - - IPhreeqcPOET(const std::string &database, const std::string &input_script, - const std::vector &solutionInitVector, - std::uint32_t n_cells) - : IPhreeqc(), n_cells(n_cells), solutionInitVector(solutionInitVector) { - this->LoadDatabaseString(database.c_str()); - this->RunString(input_script.c_str()); - - if (n_cells > 1) { - std::string copy_string = - "COPY cell 1 2-" + std::to_string(n_cells) + "\n"; - this->RunString(copy_string.c_str()); - } - } - - void queueCell(std::vector &values) { - if (this->queued_cell_pointer > this->n_cells) { - return; - } - - setCell(this->queued_cell_pointer++, values); - } - - void dequeueCells(std::vector> &values) { - if (this->queued_cell_pointer == 1) { - return; - } - - for (std::size_t i = 1; i < this->queued_cell_pointer; i++) { - std::vector cell_values; - getCell(i, cell_values); - values.push_back(cell_values); - } - - this->queued_cell_pointer = 1; - } - - void runQueuedCells(double time_step) { - if (this->queued_cell_pointer == 1) { - return; - } - run(1, this->queued_cell_pointer - 1, time_step); - } - - void setCell(int cell_number, std::vector &values) { - this->set_essential_values(cell_number, this->solutionInitVector, values); - } - - void getCell(int cell_number, std::vector &values) { - values = this->get_essential_values(cell_number, this->solutionInitVector); - } - - void run(int start_cell, int end_cell, double time_step) { - std::stringstream time_ss; - time_ss << std::fixed << std::setprecision(20) << time_step; - - const std::string runs_string = - "RUN_CELLS\n -cells " + std::to_string(start_cell) + "-" + - std::to_string(end_cell) + "\n -time_step " + time_ss.str(); - this->RunString(runs_string.c_str()); - } - - struct PhreeqcMat { - std::vector names; - std::vector ids; - std::vector> values; - }; - - PhreeqcMat getPhreeqcMat(); - - std::map raw_dumps() { - std::map dumps; - - this->SetDumpStringOn(true); - - for (const auto &[sol_id, unused] : this->raw_initials) { - std::string call_string = "DUMP\n -cells " + std::to_string(sol_id); - this->RunString(call_string.c_str()); - dumps[sol_id] = this->GetDumpString(); - } - - this->SetDumpStringOn(false); - - return dumps; - } - - using essential_names = std::array, 5>; - using ModulesArray = std::array; - - ModulesArray getModuleSizes(const std::vector &cell_ids); - - std::vector getSolutionNames(int cell_id) { - return this->raw_initials[cell_id].first[POET_SOL]; - } - -private: - // required only for simulation - essential_names initial_names; - std::uint32_t n_cells; - std::vector solutionInitVector; - - void parseInitValues(); - void parseInit(); - - essential_names dump_essential_names(std::size_t cell_number); - - cxxSolution *Get_solution(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n); - } - - cxxExchange *Get_exchange(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n); - } - - cxxKinetics *Get_kinetic(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n); - } - - cxxPPassemblage *Get_equilibrium(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(), - n); - } - - cxxSurface *Get_surface(std::size_t n) { - return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n); - } - - using RawMap = std::map>>; - - essential_names union_raws(const RawMap &raws); - - std::vector> - conc_from_essentials(const RawMap &raws, const essential_names &names); - - // void valuesFromModule(const std::string &module_name, int cell_number, - // essential_names &names, std::vector &values); - - std::string subExchangeName(std::string name) { - for (const auto &species : this->PhreeqcPtr->Get_species_list()) { - const std::string &species_name = species.s->name; - // check if `name` is a prefix of `species_name` - if (species_name.compare(0, name.size(), name) == 0) { - return species_name; - } - } - return name; - }; - - std::vector - get_essential_values(std::size_t cell_number, - const std::vector &order); - - void set_essential_values(std::size_t cell_number, - const std::vector &order, - std::vector &values); - - RawMap raw_initials; - - std::size_t queued_cell_pointer = 1; -}; \ No newline at end of file diff --git a/poet/include/KineticWrapper.hpp b/poet/include/KineticWrapper.hpp new file mode 100644 index 00000000..ba1aeda7 --- /dev/null +++ b/poet/include/KineticWrapper.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include "KineticsComp.h" +#include "WrapperBase.hpp" +#include "cxxKinetics.h" +#include +#include + +class KineticWrapper : public WrapperBase { + +public: + KineticWrapper(cxxKinetics *kinetics, + const std::vector &kin_comps); + + void get(std::span &data) const override; + void set(const std::span &data) override; + + static std::vector names(cxxKinetics *kinetics, + std::vector &kin_comps); + +private: + cxxKinetics *kinetics; + + class KineticCompWrapper : public WrapperBase { + public: + KineticCompWrapper(cxxKineticsComp &comp); + + void get(std::span &data) const override; + void set(const std::span &data) override; + + static std::vector names(const cxxKineticsComp &comp); + + private: + cxxKineticsComp &kin_comp; + }; + + std::vector> kinetic_comps; +}; \ No newline at end of file diff --git a/poet/include/POETInit.hpp b/poet/include/POETInit.hpp new file mode 100644 index 00000000..2bdc12e1 --- /dev/null +++ b/poet/include/POETInit.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include +#include + +struct POETInitCell { + std::vector solutions; + std::vector solution_primaries; + std::vector exchanger; + std::vector kinetics; + std::vector equilibrium; + std::vector surface_comps; + std::vector surface_charges; +}; + +struct POETConfig { + std::string database; + std::string input_script; + POETInitCell cell; +}; \ No newline at end of file diff --git a/poet/include/PhreeqcEngine.hpp b/poet/include/PhreeqcEngine.hpp new file mode 100644 index 00000000..3b131aa6 --- /dev/null +++ b/poet/include/PhreeqcEngine.hpp @@ -0,0 +1,96 @@ + +#pragma once + +#include "EquilibriumWrapper.hpp" +#include "ExchangeWrapper.hpp" +#include "KineticWrapper.hpp" +#include "POETInit.hpp" +#include "SolutionWrapper.hpp" +#include "SurfaceWrapper.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class PhreeqcEngine : public IPhreeqc { +public: + PhreeqcEngine(const POETConfig &config) : IPhreeqc() { + this->LoadDatabaseString(config.database.c_str()); + this->RunString(config.input_script.c_str()); + + this->init_wrappers(config.cell); + } + + void runCell(std::vector &cell_values, double time_step); + + void run(int start_cell, int end_cell) { + const std::string runs_string = "RUN_CELLS\n -cells " + + std::to_string(start_cell) + "-" + + std::to_string(end_cell) + "\nEND\n"; + this->RunString(runs_string.c_str()); + } + + void run(int start_cell, int end_cell, double time_step) { + std::stringstream time_ss; + time_ss << std::fixed << std::setprecision(20) << time_step; + + const std::string runs_string = + "RUN_CELLS\n -cells " + std::to_string(start_cell) + "-" + + std::to_string(end_cell) + "\n -time_step " + time_ss.str() + "\nEND\n"; + this->RunString(runs_string.c_str()); + } + +private: + void init_wrappers(const POETInitCell &cell); + + cxxSolution *Get_solution(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n); + } + + cxxExchange *Get_exchange(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n); + } + + cxxKinetics *Get_kinetic(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n); + } + + cxxPPassemblage *Get_equilibrium(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(), + n); + } + + cxxSurface *Get_surface(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n); + } + + void get_essential_values(std::span &data); + + void set_essential_values(const std::span &data); + + std::unique_ptr solutionWrapperPtr; + std::unique_ptr exchangeWrapperPtr; + std::unique_ptr kineticsWrapperPtr; + std::unique_ptr equilibriumWrapperPtr; + std::unique_ptr surfaceWrapperPtr; + + bool has_exchange = false; + bool has_kinetics = false; + bool has_equilibrium = false; + bool has_surface = false; +}; \ No newline at end of file diff --git a/poet/include/PhreeqcInit.hpp b/poet/include/PhreeqcInit.hpp new file mode 100644 index 00000000..b81b2d56 --- /dev/null +++ b/poet/include/PhreeqcInit.hpp @@ -0,0 +1,156 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Exchange.h" +#include "PPassemblage.h" +#include "Solution.h" +#include "Surface.h" +#include "cxxKinetics.h" + +enum { POET_SOL = 0, POET_EXCH, POET_KIN, POET_EQUIL, POET_SURF }; + +class PhreeqcInit : public IPhreeqc { +public: + PhreeqcInit(const std::string &database, const std::string &input_script); + + struct PhreeqcMat { + std::vector names; + std::vector ids; + std::vector> values; + }; + + PhreeqcMat getPhreeqcMat() const { return pqc_mat; }; + + std::map raw_dumps(); + + using essential_names = std::array, 5>; + using ModulesArray = std::array; + + ModulesArray getModuleSizes(const std::vector &cell_ids); + + std::vector getSolutionPrimaries() { + return std::vector(this->surface_primaries.begin(), + this->surface_primaries.end()); + } + + std::vector getSolutionNames(int cell_id) { + return this->raw_initials[cell_id].first[POET_SOL]; + } + + std::vector getExchanger(int id) { + if (this->exchanger.contains(id)) { + return this->exchanger[id]; + } + + return {}; + } + + std::vector getKineticsNames(int id) { + if (this->kinetics.contains(id)) { + return this->kinetics[id]; + } + + return {}; + } + + std::vector getEquilibriumNames(int id) { + if (this->equilibrium.contains(id)) { + return this->equilibrium[id]; + } + + return {}; + } + + std::vector getSurfaceCompNames(int id) { + if (this->surface_comps.contains(id)) { + return this->surface_comps[id]; + } + + return {}; + } + + std::vector getSurfaceChargeNames(int id) { + if (this->surface_charge.contains(id)) { + return this->surface_charge[id]; + } + + return {}; + } + +private: + // required only for simulation + essential_names initial_names; + + void parseInitValues(); + void parseInit(); + + std::vector dump_solution_names(int cell_number); + void dump_reactant_names(int cell_number, + const std::vector union_sol_names, + essential_names &names); + + std::vector + find_all_valence_states(const std::vector &solution_names, + const std::size_t offset); + + cxxSolution *Get_solution(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_solution_map(), n); + } + + cxxExchange *Get_exchange(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_exchange_map(), n); + } + + cxxKinetics *Get_kinetic(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_kinetics_map(), n); + } + + cxxPPassemblage *Get_equilibrium(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_pp_assemblage_map(), + n); + } + + cxxSurface *Get_surface(std::size_t n) { + return Utilities::Rxn_find(this->PhreeqcPtr->Get_Rxn_surface_map(), n); + } + + PhreeqcMat pqc_mat; + + PhreeqcMat buildPhreeqcMat(); + + std::map> exchanger; + std::map> kinetics; + std::map> equilibrium; + std::map> surface_comps; + std::map> surface_charge; + + std::set surface_primaries; + + using RawMap = std::map>>; + + essential_names union_raws(const RawMap &raws); + + std::vector> + conc_from_essentials(const RawMap &raws, const essential_names &names); + + // void valuesFromModule(const std::string &module_name, int cell_number, + // essential_names &names, std::vector &values); + + std::string subExchangeName(std::string name); + + std::vector + get_essential_values_init(std::size_t cell_number, + const std::vector &order); + + RawMap raw_initials; +}; \ No newline at end of file diff --git a/poet/include/SolutionWrapper.hpp b/poet/include/SolutionWrapper.hpp new file mode 100644 index 00000000..c75c7643 --- /dev/null +++ b/poet/include/SolutionWrapper.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include "NameDouble.h" +#include "Solution.h" +#include "WrapperBase.hpp" +#include +#include + +class SolutionWrapper : public WrapperBase { +public: + SolutionWrapper(cxxSolution *soln, + const std::vector &solution_order); + + void get(std::span &data) const; + + void set(const std::span &data); + + static std::vector + names(cxxSolution *solution, std::vector &solution_order); + +private: + cxxSolution *solution; + const std::vector solution_order; + + static constexpr std::size_t NUM_ESSENTIALS = 3; +}; \ No newline at end of file diff --git a/poet/include/SurfaceWrapper.hpp b/poet/include/SurfaceWrapper.hpp new file mode 100644 index 00000000..b96b902a --- /dev/null +++ b/poet/include/SurfaceWrapper.hpp @@ -0,0 +1,71 @@ +#pragma once + +#include "Surface.h" +#include "SurfaceCharge.h" +#include "SurfaceComp.h" +#include "WrapperBase.hpp" +#include "phrqtype.h" +#include +#include +#include +#include +#include +#include + +class SurfaceWrapper : public WrapperBase { +public: + SurfaceWrapper(cxxSurface *surf, + const std::set &solution_primaries, + const std::vector &comp_formulas, + const std::vector &charge_names); + + void get(std::span &surface) const override; + + void set(const std::span &surface) override; + + static std::vector + names(cxxSurface *surface, const std::set &solution_primaries, + std::vector &comp_formulas, + std::vector &charge_names); + +private: + cxxSurface *surface; + + class SurfaceCompWrapper : public WrapperBase { + public: + SurfaceCompWrapper(cxxSurfaceComp &comp); + + void get(std::span &surface) const; + + void set(const std::span &surface); + + static std::vector names(cxxSurfaceComp &comp); + + private: + cxxSurfaceComp &surface_comp; + static constexpr std::size_t NUM_NOT_TOTALS = 3; + + std::vector total_names; + }; + + class SurfaceChargeWrapper : public WrapperBase { + public: + SurfaceChargeWrapper(cxxSurfaceCharge &charge, + const std::set &solution_primaries); + void get(std::span &surface) const; + + void set(const std::span &surface); + + static std::vector + names(cxxSurfaceCharge *charge, const std::set &primaries); + + private: + cxxSurfaceCharge &surface_charge; + static constexpr std::size_t NUM_NOT_TOTALS = 5; + + std::set primaries; + }; + + std::vector> surface_comps; + std::vector> surface_charges; +}; \ No newline at end of file diff --git a/poet/include/WrapperBase.hpp b/poet/include/WrapperBase.hpp new file mode 100644 index 00000000..b12d36eb --- /dev/null +++ b/poet/include/WrapperBase.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include +#include + +class WrapperBase { +public: + virtual ~WrapperBase() = default; + + std::size_t size() const { return this->num_elements; }; + + virtual void get(std::span &data) const = 0; + + virtual void set(const std::span &data) = 0; + +protected: + std::size_t num_elements = 0; +}; \ No newline at end of file diff --git a/poet/src/Engine.cpp b/poet/src/Engine.cpp new file mode 100644 index 00000000..2f327bbf --- /dev/null +++ b/poet/src/Engine.cpp @@ -0,0 +1,220 @@ +#include "PhreeqcEngine.hpp" +#include +#include + +void PhreeqcEngine::init_wrappers(const POETInitCell &cell) { + + // TODO: Implement the rest of the wrappers. Currently supported: EXCHANGE + + // Solutions + this->solutionWrapperPtr = + std::make_unique(this->Get_solution(1), cell.solutions); + + if (this->Get_exchange(1) != nullptr) { + this->exchangeWrapperPtr = std::make_unique( + this->Get_exchange(1), cell.exchanger); + this->has_exchange = true; + } + + if (this->Get_kinetic(1) != nullptr) { + this->kineticsWrapperPtr = + std::make_unique(this->Get_kinetic(1), cell.kinetics); + + this->has_kinetics = true; + } + + if (this->Get_equilibrium(1) != nullptr) { + this->equilibriumWrapperPtr = std::make_unique( + this->Get_equilibrium(1), cell.equilibrium); + + this->has_equilibrium = true; + } + + if (this->Get_surface(1) != nullptr) { + std::set primaries(cell.solution_primaries.begin(), + cell.solution_primaries.end()); + this->surfaceWrapperPtr = std::make_unique( + this->Get_surface(1), primaries, cell.surface_comps, + cell.surface_charges); + + this->has_surface = true; + } +} + +void PhreeqcEngine::get_essential_values(std::span &data) { + + this->solutionWrapperPtr->get(data); + + std::size_t offset = this->solutionWrapperPtr->size(); + + if (this->has_exchange) { + std::span exch_span{ + data.subspan(offset, this->exchangeWrapperPtr->size())}; + this->exchangeWrapperPtr->get(exch_span); + + offset += this->exchangeWrapperPtr->size(); + } + + if (this->has_kinetics) { + std::span kin_span{ + data.subspan(offset, this->kineticsWrapperPtr->size())}; + this->kineticsWrapperPtr->get(kin_span); + + offset += this->kineticsWrapperPtr->size(); + } + + if (this->has_equilibrium) { + std::span equ_span{ + data.subspan(offset, this->equilibriumWrapperPtr->size())}; + this->equilibriumWrapperPtr->get(equ_span); + + offset += this->equilibriumWrapperPtr->size(); + } + + if (this->has_surface) { + std::span surf_span{ + data.subspan(offset, this->surfaceWrapperPtr->size())}; + this->surfaceWrapperPtr->get(surf_span); + } + + // // Exchange + // if (this->Get_exchange(cell_number) != NULL) { + // std::vector exch_values(this->exchangeWrapperPtr->size()); + + // std::span exch_span{exch_values}; + + // this->exchangeWrapperPtr->get(exch_span); + + // // this->Get_exchange(cell_number)->get_essential_values(exch_values); + // essentials.insert(essentials.end(), exch_values.begin(), + // exch_values.end()); + // } + + // // Kinetics + // if (this->Get_kinetic(cell_number) != NULL) { + // std::vector kin_values; + + // this->Get_kinetic(cell_number)->get_essential_values(kin_values); + // essentials.insert(essentials.end(), kin_values.begin(), + // kin_values.end()); + // } + + // // PPassemblage + // if (this->Get_equilibrium(cell_number) != NULL) { + // std::vector equ_values; + + // this->Get_equilibrium(cell_number)->get_essential_values(equ_values); + // essentials.insert(essentials.end(), equ_values.begin(), + // equ_values.end()); + // } + + // // Surface + // if (this->Get_surface(cell_number) != NULL) { + // std::vector surf_values(this->surfaceWrapperPtr->size()); + // std::span surf_span{surf_values}; + + // this->surfaceWrapperPtr->get(surf_span); + + // // this->Get_surface(cell_number)->get_essential_values(surf_values); + // essentials.insert(essentials.end(), surf_values.begin(), + // surf_values.end()); + // } + + // return essentials; +} + +void PhreeqcEngine::set_essential_values(const std::span &data) { + + this->solutionWrapperPtr->set(data); + this->PhreeqcPtr->initial_solutions_poet(1); + + std::size_t offset = this->solutionWrapperPtr->size(); + + if (this->has_exchange) { + std::span exch_span{ + data.subspan(offset, this->exchangeWrapperPtr->size())}; + this->exchangeWrapperPtr->set(exch_span); + + offset += this->exchangeWrapperPtr->size(); + } + + if (this->has_kinetics) { + std::span kin_span{ + data.subspan(offset, this->kineticsWrapperPtr->size())}; + this->kineticsWrapperPtr->set(kin_span); + + offset += this->kineticsWrapperPtr->size(); + } + + if (this->has_equilibrium) { + std::span equ_span{ + data.subspan(offset, this->equilibriumWrapperPtr->size())}; + this->equilibriumWrapperPtr->set(equ_span); + + offset += this->equilibriumWrapperPtr->size(); + } + + if (this->has_surface) { + std::span surf_span{ + data.subspan(offset, this->surfaceWrapperPtr->size())}; + this->surfaceWrapperPtr->set(surf_span); + } + + // // Solutions + // if (this->Get_solution(cell_number) != NULL) { + // this->Get_solution(cell_number)->set_essential_values(dump_it, order); + // this->PhreeqcPtr->initial_solutions_poet(cell_number); + // } + + // // Exchange + // if (this->Get_exchange(cell_number) != NULL) { + // auto &exch_pointer = this->exchangeWrapperPtr; + // std::span exch_span{dump_it, exch_pointer->size()}; + // exch_pointer->set(exch_span); + + // dump_it += exch_pointer->size(); + + // // this->Get_exchange(cell_number)->set_essential_values(dump_it); + // // this->PhreeqcPtr->Rxn_new_exchange.insert(cell_number); + // // this->Get_exchange(cell_number)->Set_new_def(true); + // // this->Get_exchange(cell_number)->Set_solution_equilibria(true); + // // this->Get_exchange(cell_number)->Set_n_solution(cell_number); + // // this->PhreeqcPtr->initial_exchangers(FALSE); + // } + + // // Kinetics + // if (this->Get_kinetic(cell_number) != NULL) { + // this->Get_kinetic(cell_number)->set_essential_values(dump_it); + // } + + // // PPassemblage + // if (this->Get_equilibrium(cell_number) != NULL) { + // this->Get_equilibrium(cell_number)->set_essential_values(dump_it); + // } + + // // Surface + // if (this->Get_surface(cell_number) != NULL) { + // auto *tmp = this->Get_surface(cell_number); + // auto &surf_pointer = this->surfaceWrapperPtr; + // std::span surf_span{dump_it, surf_pointer->size()}; + // surf_pointer->set(surf_span); + + // dump_it += surf_pointer->size(); + + // // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number); + + // // this->Get_surface(cell_number)->set_essential_values(dump_it); + // // this->PhreeqcPtr->Rxn_new_surface.insert(cell_number); + // // this->PhreeqcPtr->initial_surfaces(TRUE); + // } +} + +void PhreeqcEngine::runCell(std::vector &cell_values, + double time_step) { + // skip ID + std::span cell_data{cell_values.begin() + 1, cell_values.end()}; + + this->set_essential_values(cell_data); + this->run(1, 1, time_step); + this->get_essential_values(cell_data); +} \ No newline at end of file diff --git a/poet/src/EquilibriumCompWrapper.cpp b/poet/src/EquilibriumCompWrapper.cpp new file mode 100644 index 00000000..2a7c5672 --- /dev/null +++ b/poet/src/EquilibriumCompWrapper.cpp @@ -0,0 +1,30 @@ +#include "EquilibriumWrapper.hpp" + +EquilibriumWrapper::EquilibriumCompWrapper::EquilibriumCompWrapper( + cxxPPassemblageComp &comp_) + : comp(comp_) { + this->num_elements = 2; +} + +void EquilibriumWrapper::EquilibriumCompWrapper::get( + std::span &data) const { + data[0] = this->comp.Get_moles(); + data[1] = this->comp.Get_si(); +} + +void EquilibriumWrapper::EquilibriumCompWrapper::set( + const std::span &data) { + this->comp.Set_moles(data[0]); + this->comp.Set_si(data[1]); +} + +std::vector EquilibriumWrapper::EquilibriumCompWrapper::names( + const cxxPPassemblageComp &comp) { + std::vector names; + + const std::string &comp_name = comp.Get_name(); + names.push_back(comp_name); + names.push_back(comp_name + "_si"); + + return names; +} \ No newline at end of file diff --git a/poet/src/EquilibriumWrapper.cpp b/poet/src/EquilibriumWrapper.cpp new file mode 100644 index 00000000..e5513569 --- /dev/null +++ b/poet/src/EquilibriumWrapper.cpp @@ -0,0 +1,66 @@ +#include "EquilibriumWrapper.hpp" + +#include +#include + +EquilibriumWrapper::EquilibriumWrapper( + cxxPPassemblage *ppassemblage, + const std::vector &ppassemblage_names) + : ppassemblage(ppassemblage) { + for (const auto &ppassemblage_name : ppassemblage_names) { + auto it = std::find_if( + ppassemblage->Get_pp_assemblage_comps().begin(), + ppassemblage->Get_pp_assemblage_comps().end(), + [&](const auto &comp) { return comp.first == ppassemblage_name; }); + + if (it == ppassemblage->Get_pp_assemblage_comps().end()) { + throw std::runtime_error( + "Equilibrium component not found in Phreeqc variables"); + } + + equilibrium_comps.push_back( + std::make_unique(it->second)); + + num_elements += equilibrium_comps.back()->size(); + } +} + +void EquilibriumWrapper::get(std::span &data) const { + std::size_t offset = 0; + + for (const auto &comp : equilibrium_comps) { + std::span comp_span = data.subspan(offset, comp->size()); + + comp->get(comp_span); + + offset += comp->size(); + } +} + +void EquilibriumWrapper::set(const std::span &data) { + std::size_t offset = 0; + + for (const auto &comp : equilibrium_comps) { + std::span comp_span = data.subspan(offset, comp->size()); + + comp->set(comp_span); + + offset += comp->size(); + } +} + +std::vector +EquilibriumWrapper::names(const cxxPPassemblage *ppassemblage, + std::vector &ppassemblage_names) { + std::vector names; + + for (const auto &comp : ppassemblage->Get_pp_assemblage_comps()) { + std::vector comp_names = + EquilibriumCompWrapper::names(comp.second); + names.insert(names.end(), comp_names.begin(), comp_names.end()); + + ppassemblage_names.push_back(comp.first); + } + + return names; +} \ No newline at end of file diff --git a/poet/src/ExchangeCompWrapper.cpp b/poet/src/ExchangeCompWrapper.cpp new file mode 100644 index 00000000..a259b0e3 --- /dev/null +++ b/poet/src/ExchangeCompWrapper.cpp @@ -0,0 +1,41 @@ +#include "ExchangeWrapper.hpp" + +ExchangeWrapper::ExchangeCompWrapper::ExchangeCompWrapper(cxxExchComp &comp) + : exch_comp(comp) { + const std::size_t totals_size = exch_comp.Get_totals().size() - 1; + this->num_elements = totals_size + NUM_NOT_TOTALS; +} + +void ExchangeWrapper::ExchangeCompWrapper::get( + std::span &exchange) const { + exchange[0] = exch_comp.Get_totals().find(exch_comp.Get_formula())->second; + exchange[1] = exch_comp.Get_charge_balance(); + exchange[2] = exch_comp.Get_la(); + exchange[3] = exch_comp.Get_phase_proportion(); + exchange[4] = exch_comp.Get_formula_z(); + + for (const auto &[name, value] : exch_comp.Get_totals()) { + if (name != exch_comp.Get_formula()) { + exchange[NUM_NOT_TOTALS + + std::distance(exch_comp.Get_totals().begin(), + exch_comp.Get_totals().find(name))] = value; + } + } +} + +void ExchangeWrapper::ExchangeCompWrapper::set( + const std::span &exchange) { + exch_comp.Get_totals().find(exch_comp.Get_formula())->second = exchange[0]; + exch_comp.Set_charge_balance(exchange[1]); + exch_comp.Set_la(exchange[2]); + exch_comp.Set_phase_proportion(exchange[3]); + exch_comp.Set_formula_z(exchange[4]); + + for (auto &[name, value] : exch_comp.Get_totals()) { + if (name != exch_comp.Get_formula()) { + value = exchange[NUM_NOT_TOTALS + + std::distance(exch_comp.Get_totals().begin(), + exch_comp.Get_totals().find(name))]; + } + } +} \ No newline at end of file diff --git a/poet/src/ExchangeWrapper.cpp b/poet/src/ExchangeWrapper.cpp new file mode 100644 index 00000000..378e10e0 --- /dev/null +++ b/poet/src/ExchangeWrapper.cpp @@ -0,0 +1,101 @@ +#include "ExchangeWrapper.hpp" + +ExchangeWrapper::ExchangeWrapper(cxxExchange *exch, + const std::vector &exchanger) + : exchange(exch) { + + for (const auto &comp_name : exchanger) { + auto it = std::find_if(exchange->Get_exchange_comps().begin(), + exchange->Get_exchange_comps().end(), + [&](const cxxExchComp &comp) { + return comp.Get_formula() == comp_name; + }); + + if (it == exchange->Get_exchange_comps().end()) { + throw std::runtime_error( + "Exchange component not found in Phreeqc variables"); + } + + exchange_comps.push_back(std::make_unique(*it)); + + num_elements += exchange_comps.back()->size(); + } + + // const std::size_t defined_comps = exchange->Get_exchange_comps().size(); + + // auto header_it = remaining_field_header.begin(); + + // while (header_it != remaining_field_header.end() && + // exchange_comps.size() < defined_comps) { + // const std::string formular = *header_it; + + // auto it = std::find_if(exchange->Get_exchange_comps().begin(), + // exchange->Get_exchange_comps().end(), + // [&](const cxxExchComp &comp) { + // return comp.Get_formula() == formular; + // }); + + // if (it != exchange->Get_exchange_comps().end()) { + // const size_t i = this->exchange_comps.size(); + + // exchange_comps.push_back(std::make_unique(*it)); + // header_it += this->exchange_comps[i]->size(); + // num_elements += this->exchange_comps[i]->size(); + // continue; + // } + + // header_it++; + // } + + // if (exchange_comps.size() != defined_comps) { + // throw std::runtime_error( + // "Not all exchange components found in Phreeqc variables"); + // } +} + +void ExchangeWrapper::get(std::span &exchange) const { + std::size_t offset = 0; + + for (const auto &comp : exchange_comps) { + std::span comp_span = exchange.subspan(offset, comp->size()); + comp->get(comp_span); + offset += comp->size(); + } +} + +void ExchangeWrapper::set(const std::span &exchange) { + std::size_t offset = 0; + + for (const auto &comp : exchange_comps) { + std::span comp_span = exchange.subspan(offset, comp->size()); + comp->set(comp_span); + offset += comp->size(); + } +} + +std::vector +ExchangeWrapper::names(cxxExchange *exchange, + std::vector &exchange_formulas) { + std::vector e_names; + + for (const auto &comp : exchange->Get_exchange_comps()) { + + const std::string &formular = comp.Get_formula(); + exchange_formulas.push_back(formular); + e_names.push_back(formular); + + e_names.push_back(formular + "_cb"); + e_names.push_back(formular + "_la"); + e_names.push_back(formular + "_phase_proportion"); + e_names.push_back(formular + "_formular_z"); + + for (const auto &total : comp.Get_totals()) { + if (total.first == formular) { + continue; + } + e_names.push_back(total.first + formular); + } + } + + return e_names; +} \ No newline at end of file diff --git a/poet/src/Init.cpp b/poet/src/Extraction.cpp similarity index 63% rename from poet/src/Init.cpp rename to poet/src/Extraction.cpp index 2aaa9341..d8b1ea1f 100644 --- a/poet/src/Init.cpp +++ b/poet/src/Extraction.cpp @@ -1,4 +1,6 @@ -#include +#include + +#include "Solution.h" #include #include @@ -37,9 +39,19 @@ createConcVector(const std::vector &conc_names, return conc_vec; } -void IPhreeqcPOET::parseInit() { +PhreeqcInit::PhreeqcInit(const std::string &database, + const std::string &input_script) + : IPhreeqc() { + this->LoadDatabaseString(database.c_str()); + this->RunString(input_script.c_str()); - std::map>> init_values; + this->parseInit(); + pqc_mat = this->buildPhreeqcMat(); +} + +void PhreeqcInit::parseInit() { + + std::map init_values; for (const auto &[id, val] : this->PhreeqcPtr->Get_Rxn_solution_map()) { // A key less than zero indicates an internal solution @@ -47,16 +59,50 @@ void IPhreeqcPOET::parseInit() { continue; } - const essential_names curr_conc_names = this->dump_essential_names(id); + essential_names curr_conc_names; + + curr_conc_names[POET_SOL] = dump_solution_names(id); + init_values[id] = curr_conc_names; + + // auto pair = + // std::make_pair(curr_conc_names, + // this->get_essential_values_init(id, + // curr_conc_names[0])); + // raw_initials[id] = pair; + } + + std::vector union_solution_names; + + for (const auto &[id, val] : init_values) { + union_solution_names = + unionStringVectors(union_solution_names, val[POET_SOL]); + } + + for (auto &[id, val] : init_values) { + dump_reactant_names(id, union_solution_names, val); + + std::vector solution_order(val[POET_SOL].begin() + 3, + val[POET_SOL].end()); auto pair = std::make_pair( - curr_conc_names, this->get_essential_values(id, curr_conc_names[0])); + val, this->get_essential_values_init(id, solution_order)); raw_initials[id] = pair; } } -IPhreeqcPOET::essential_names IPhreeqcPOET::union_raws(const RawMap &raws) { - IPhreeqcPOET::essential_names names; +std::string PhreeqcInit::subExchangeName(std::string name) { + for (const auto &species : this->PhreeqcPtr->Get_species_list()) { + const std::string &species_name = species.s->name; + // check if `name` is a prefix of `species_name` + if (species_name.compare(0, name.size(), name) == 0) { + return species_name; + } + } + return name; +}; + +PhreeqcInit::essential_names PhreeqcInit::union_raws(const RawMap &raws) { + PhreeqcInit::essential_names names; for (const auto &[sol_id, val] : raws) { for (std::size_t i = 0; i < names.size(); i++) { names[i] = unionStringVectors(names[i], val.first[i]); @@ -71,8 +117,8 @@ IPhreeqcPOET::essential_names IPhreeqcPOET::union_raws(const RawMap &raws) { } std::vector> -IPhreeqcPOET::conc_from_essentials(const RawMap &raws, - const essential_names &names) { +PhreeqcInit::conc_from_essentials(const RawMap &raws, + const essential_names &names) { std::vector> values; for (const auto &[key, val] : raws) { @@ -94,10 +140,10 @@ IPhreeqcPOET::conc_from_essentials(const RawMap &raws, return values; } -IPhreeqcPOET::PhreeqcMat IPhreeqcPOET::getPhreeqcMat() { +PhreeqcInit::PhreeqcMat PhreeqcInit::buildPhreeqcMat() { std::vector> values; - const IPhreeqcPOET::essential_names ess_names = + const PhreeqcInit::essential_names ess_names = this->union_raws(this->raw_initials); const std::vector> conc_values = @@ -117,8 +163,8 @@ IPhreeqcPOET::PhreeqcMat IPhreeqcPOET::getPhreeqcMat() { return {.names = conc_names, .ids = ids, .values = conc_values}; } -IPhreeqcPOET::ModulesArray -IPhreeqcPOET::getModuleSizes(const std::vector &cell_ids) { +PhreeqcInit::ModulesArray +PhreeqcInit::getModuleSizes(const std::vector &cell_ids) { ModulesArray module_sizes; RawMap raws; @@ -126,7 +172,7 @@ IPhreeqcPOET::getModuleSizes(const std::vector &cell_ids) { raws[id] = this->raw_initials[id]; } - const IPhreeqcPOET::essential_names ess_names = this->union_raws(raws); + const PhreeqcInit::essential_names ess_names = this->union_raws(raws); std::transform(ess_names.begin(), ess_names.end(), module_sizes.begin(), [](const auto &mod) { return mod.size(); }); diff --git a/poet/src/GetSet.cpp b/poet/src/GetSet.cpp deleted file mode 100644 index 8a7c4aff..00000000 --- a/poet/src/GetSet.cpp +++ /dev/null @@ -1,119 +0,0 @@ -#include - -IPhreeqcPOET::essential_names -IPhreeqcPOET::dump_essential_names(std::size_t cell_number) { - - essential_names eNames; - - // Solutions - if (this->Get_solution(cell_number) != NULL) { - std::vector &eSolNames = eNames[POET_SOL]; - this->Get_solution(cell_number)->dump_essential_names(eSolNames); - } - - // Exchange - if (this->Get_exchange(cell_number) != NULL) { - std::vector &eExchNames = eNames[POET_EXCH]; - this->Get_exchange(cell_number)->dump_essential_names(eExchNames); - } - - // Kinetics - if (this->Get_kinetic(cell_number) != NULL) { - std::vector &eKinNames = eNames[POET_KIN]; - this->Get_kinetic(cell_number)->dump_essential_names(eKinNames); - } - - // PPassemblage - if (this->Get_equilibrium(cell_number) != NULL) { - std::vector &eEquNames = eNames[POET_EQUIL]; - this->Get_equilibrium(cell_number)->dump_essential_names(eEquNames); - } - - // Surface - if (this->Get_surface(cell_number) != NULL) { - std::vector &eSurfNames = eNames[POET_SURF]; - this->Get_surface(cell_number)->dump_essential_names(eSurfNames); - } - - return eNames; -} - -std::vector -IPhreeqcPOET::get_essential_values(std::size_t cell_number, - const std::vector &order) { - std::vector essentials; - - // Solutions - if (this->Get_solution(cell_number) != NULL) { - std::vector sol_values; - - this->Get_solution(cell_number)->get_essential_values(sol_values, order); - essentials.insert(essentials.end(), sol_values.begin(), sol_values.end()); - } - - // Exchange - if (this->Get_exchange(cell_number) != NULL) { - std::vector exch_values; - - this->Get_exchange(cell_number)->get_essential_values(exch_values); - essentials.insert(essentials.end(), exch_values.begin(), exch_values.end()); - } - - // Kinetics - if (this->Get_kinetic(cell_number) != NULL) { - std::vector kin_values; - - this->Get_kinetic(cell_number)->get_essential_values(kin_values); - essentials.insert(essentials.end(), kin_values.begin(), kin_values.end()); - } - - // PPassemblage - if (this->Get_equilibrium(cell_number) != NULL) { - std::vector equ_values; - - this->Get_equilibrium(cell_number)->get_essential_values(equ_values); - essentials.insert(essentials.end(), equ_values.begin(), equ_values.end()); - } - - // Surface - if (this->Get_surface(cell_number) != NULL) { - std::vector surf_values; - - this->Get_surface(cell_number)->get_essential_values(surf_values); - essentials.insert(essentials.end(), surf_values.begin(), surf_values.end()); - } - - return essentials; -} - -void IPhreeqcPOET::set_essential_values(std::size_t cell_number, - const std::vector &order, - std::vector &values) { - - auto dump_it = values.begin(); - - // Solutions - if (this->Get_solution(cell_number) != NULL) { - this->Get_solution(cell_number)->set_essential_values(dump_it, order); - } - - // Exchange - if (this->Get_exchange(cell_number) != NULL) { - this->Get_exchange(cell_number)->set_essential_values(dump_it); - } - - // Kinetics - if (this->Get_kinetic(cell_number) != NULL) { - this->Get_kinetic(cell_number)->set_essential_values(dump_it); - } - - // PPassemblage - if (this->Get_equilibrium(cell_number) != NULL) { - this->Get_equilibrium(cell_number)->set_essential_values(dump_it); - } - - // Surface - if (this->Get_surface(cell_number) != NULL) { - this->Get_surface(cell_number)->set_essential_values(dump_it); - } -} \ No newline at end of file diff --git a/poet/src/InitGetSet.cpp b/poet/src/InitGetSet.cpp new file mode 100644 index 00000000..11549f2b --- /dev/null +++ b/poet/src/InitGetSet.cpp @@ -0,0 +1,251 @@ +#include "EquilibriumWrapper.hpp" +#include "ExchangeWrapper.hpp" +#include "KineticWrapper.hpp" +#include "PhreeqcInit.hpp" +#include "SolutionWrapper.hpp" +#include "SurfaceWrapper.hpp" +#include "global_structures.h" + +#include +#include +#include +#include +#include +#include + +std::vector PhreeqcInit::dump_solution_names(int cell_number) { + constexpr std::size_t SKIP_H_O_CB = 3; + + std::vector placeholder; + + return find_all_valence_states( + SolutionWrapper::names(this->Get_solution(cell_number), placeholder), + SKIP_H_O_CB); +} + +void PhreeqcInit::dump_reactant_names( + int cell_number, const std::vector union_sol_names, + essential_names &names) { + // Exchange + if (this->Get_exchange(cell_number) != NULL) { + auto *exchange = this->Get_exchange(cell_number); + + std::vector exch_formulas; + names[POET_EXCH] = ExchangeWrapper::names(exchange, exch_formulas); + + this->exchanger[cell_number] = exch_formulas; + } + + // Kinetics + if (this->Get_kinetic(cell_number) != NULL) { + auto *kinetic = this->Get_kinetic(cell_number); + + std::vector kin_comps; + + names[POET_KIN] = KineticWrapper::names(kinetic, kin_comps); + + this->kinetics[cell_number] = kin_comps; + } + + // PPassemblage + if (this->Get_equilibrium(cell_number) != NULL) { + auto *equilibrium = this->Get_equilibrium(cell_number); + + std::vector equ_names; + + names[POET_EQUIL] = EquilibriumWrapper::names(equilibrium, equ_names); + + this->equilibrium[cell_number] = equ_names; + } + + // Surface + if (this->Get_surface(cell_number) != NULL) { + auto *surface = this->Get_surface(cell_number); + + if (this->surface_primaries.empty()) { + for (std::size_t i = 3; i < union_sol_names.size(); i++) { + const auto master_primary = this->PhreeqcPtr->master_bsearch_primary( + union_sol_names[i].c_str()); + if (master_primary != NULL) { + this->surface_primaries.insert(master_primary->elt->name); + } + } + } + + std::vector comp_formulas; + std::vector charge_names; + + names[POET_SURF] = SurfaceWrapper::names(surface, this->surface_primaries, + comp_formulas, charge_names); + + this->surface_comps[cell_number] = comp_formulas; + this->surface_charge[cell_number] = charge_names; + } +} + +std::vector PhreeqcInit::find_all_valence_states( + const std::vector &solution_names, const std::size_t offset) { + std::vector solution_with_valences( + solution_names.begin(), solution_names.begin() + offset); + + // to avoid duplicates store already evaluated master species + std::set master_species_found; + + for (std::size_t i = offset; i < solution_names.size(); i++) { + const auto master_primary = + this->PhreeqcPtr->master_bsearch_primary(solution_names[i].c_str()); + + assert(master_primary != NULL); + + const std::string master_species(master_primary->elt->name); + + // in case we already found valences for this master species we skip it + if (master_species_found.contains(master_species)) { + continue; + } + + master_species_found.insert(master_species); + + bool has_valences = false; + std::size_t inner_loop_j = 0; + std::size_t last_valence = inner_loop_j; + + // we HAVE to assume master species are already sorted! + // loop over all known master species + for (inner_loop_j = 0; inner_loop_j < this->PhreeqcPtr->master.size(); + inner_loop_j++) { + std::string curr_master_species( + this->PhreeqcPtr->master[inner_loop_j]->elt->name); + + if (curr_master_species == master_species) { + last_valence = inner_loop_j; + + while (last_valence < this->PhreeqcPtr->master.size() - 1) { + std::string next_master_species( + this->PhreeqcPtr->master[last_valence + 1]->elt->name); + + // check if the next name starts with the current master species + if (next_master_species.compare(0, master_species.size(), + master_species) != 0) { + break; + } + + // check if the next character is an opening parenthesis + if (next_master_species[master_species.size()] != '(') { + break; + } + + // if this is all true we have a valence state + has_valences = true; + last_valence++; + } + break; + } + } + + // in case we found valences we add them to the solution + if (has_valences) { + // skip the master species and only add the valences + for (std::size_t k = inner_loop_j + 1; k <= last_valence; k++) { + solution_with_valences.push_back( + this->PhreeqcPtr->master[k]->elt->name); + } + continue; + } + + // otherwise we just add the master species without any valences + solution_with_valences.push_back(master_species); + } + + return solution_with_valences; +} + +std::vector +PhreeqcInit::get_essential_values_init(std::size_t cell_number, + const std::vector &order) { + std::vector essentials; + + // Solutions + if (this->Get_solution(cell_number) != NULL) { + SolutionWrapper solution(this->Get_solution(cell_number), order); + std::vector sol_values(solution.size()); + std::span sol_span{sol_values}; + + solution.get(sol_span); + + essentials.insert(essentials.end(), sol_values.begin(), sol_values.end()); + } + + // Exchange + if (this->Get_exchange(cell_number) != NULL) { + ExchangeWrapper exchange(this->Get_exchange(cell_number), + this->exchanger[cell_number]); + + std::vector exch_values(exchange.size()); + std::span exch_span{exch_values}; + + exchange.get(exch_span); + + essentials.insert(essentials.end(), exch_values.begin(), exch_values.end()); + } + + // Kinetics + if (this->Get_kinetic(cell_number) != NULL) { + KineticWrapper kinetic(this->Get_kinetic(cell_number), + this->kinetics[cell_number]); + + std::vector kin_values(kinetic.size()); + std::span kin_span{kin_values}; + + kinetic.get(kin_span); + + essentials.insert(essentials.end(), kin_values.begin(), kin_values.end()); + } + + // PPassemblage + if (this->Get_equilibrium(cell_number) != NULL) { + EquilibriumWrapper equilibrium(this->Get_equilibrium(cell_number), + this->equilibrium[cell_number]); + + std::vector equ_values(equilibrium.size()); + std::span equ_span{equ_values}; + + equilibrium.get(equ_span); + + essentials.insert(essentials.end(), equ_values.begin(), equ_values.end()); + } + + // Surface + if (this->Get_surface(cell_number) != NULL) { + + SurfaceWrapper surface( + this->Get_surface(cell_number), this->surface_primaries, + this->surface_comps[cell_number], this->surface_charge[cell_number]); + + std::vector surf_values(surface.size()); + std::span surf_span{surf_values}; + + surface.get(surf_span); + + essentials.insert(essentials.end(), surf_values.begin(), surf_values.end()); + } + + return essentials; +} + +std::map PhreeqcInit::raw_dumps() { + std::map dumps; + + this->SetDumpStringOn(true); + + for (const auto &[sol_id, unused] : this->raw_initials) { + std::string call_string = + "DUMP\n -cells " + std::to_string(sol_id) + "\nEND\n"; + this->RunString(call_string.c_str()); + dumps[sol_id] = this->GetDumpString(); + } + + this->SetDumpStringOn(false); + + return dumps; +} \ No newline at end of file diff --git a/poet/src/KineticCompWrapper.cpp b/poet/src/KineticCompWrapper.cpp new file mode 100644 index 00000000..ebfec6e6 --- /dev/null +++ b/poet/src/KineticCompWrapper.cpp @@ -0,0 +1,40 @@ +#include "KineticWrapper.hpp" +#include +#include +#include + +KineticWrapper::KineticCompWrapper::KineticCompWrapper(cxxKineticsComp &comp) + : kin_comp(comp) { + this->num_elements = kin_comp.Get_d_params().size() + 1; +} + +void KineticWrapper::KineticCompWrapper::get(std::span &data) const { + data[0] = this->kin_comp.Get_m(); + + std::size_t next_index = 1; + for (const auto ¶m : this->kin_comp.Get_d_params()) { + data[next_index++] = param; + } +} +void KineticWrapper::KineticCompWrapper::set(const std::span &data) { + this->kin_comp.Set_m(data[0]); + + std::size_t next_index = 1; + for (auto ¶m : this->kin_comp.Get_d_params()) { + param = data[next_index++]; + } +} + +std::vector +KineticWrapper::KineticCompWrapper::names(const cxxKineticsComp &comp) { + std::vector names; + + const std::string &comp_name = comp.Get_rate_name(); + names.push_back(comp_name); + + for (std::size_t i = 0; i < comp.Get_d_params().size(); i++) { + names.push_back(comp_name + "_p" + std::to_string(i + 1)); + } + + return names; +} \ No newline at end of file diff --git a/poet/src/KineticWrapper.cpp b/poet/src/KineticWrapper.cpp new file mode 100644 index 00000000..9a2f85e3 --- /dev/null +++ b/poet/src/KineticWrapper.cpp @@ -0,0 +1,78 @@ +#include "KineticWrapper.hpp" +#include + +KineticWrapper::KineticWrapper(cxxKinetics *kinetics, + const std::vector &kin_comps) + : kinetics(kinetics) { + + for (const auto &kin_name : kin_comps) { + auto it = std::find_if(kinetics->Get_kinetics_comps().begin(), + kinetics->Get_kinetics_comps().end(), + [&](const cxxKineticsComp &comp) { + return comp.Get_rate_name() == kin_name; + }); + + if (it == kinetics->Get_kinetics_comps().end()) { + throw std::runtime_error( + "Kinetic component not found in Phreeqc variables"); + } + + kinetic_comps.push_back(std::make_unique(*it)); + + num_elements += kinetic_comps.back()->size(); + } +} + +void KineticWrapper::get(std::span &data) const { + std::size_t offset = 0; + + for (const auto &comp : kinetic_comps) { + std::span comp_span = data.subspan(offset, comp->size()); + + comp->get(comp_span); + + offset += comp->size(); + } +} + +void KineticWrapper::set(const std::span &data) { + std::size_t offset = 0; + + for (const auto &comp : kinetic_comps) { + std::span comp_span = data.subspan(offset, comp->size()); + + comp->set(comp_span); + + offset += comp->size(); + } +} + +std::vector +KineticWrapper::names(cxxKinetics *kinetics, + std::vector &kin_comps) { + std::vector names; + + for (const auto &comp : kinetics->Get_kinetics_comps()) { + std::vector comp_names = KineticCompWrapper::names(comp); + names.insert(names.end(), comp_names.begin(), comp_names.end()); + + kin_comps.push_back(comp.Get_rate_name()); + } + + for (const auto &kin_name : kin_comps) { + auto it = std::find_if(kinetics->Get_kinetics_comps().begin(), + kinetics->Get_kinetics_comps().end(), + [&](const cxxKineticsComp &comp) { + return comp.Get_rate_name() == kin_name; + }); + + if (it == kinetics->Get_kinetics_comps().end()) { + throw std::runtime_error( + "Kinetic component not found in Phreeqc variables"); + } + + names.push_back(it->Get_rate_name()); + } + + return names; +} \ No newline at end of file diff --git a/poet/src/SolutionWrapper.cpp b/poet/src/SolutionWrapper.cpp new file mode 100644 index 00000000..eb22a175 --- /dev/null +++ b/poet/src/SolutionWrapper.cpp @@ -0,0 +1,63 @@ +#include "SolutionWrapper.hpp" +#include "NameDouble.h" +#include + +SolutionWrapper::SolutionWrapper( + cxxSolution *soln, const std::vector &solution_order_) + : solution(soln), solution_order(solution_order_) { + this->num_elements = solution_order.size() + NUM_ESSENTIALS; + + auto &totals = solution->Get_totals(); +} + +void SolutionWrapper::get(std::span &data) const { + data[0] = solution->Get_total_h(); + data[1] = solution->Get_total_o(); + data[2] = solution->Get_cb(); + + std::size_t i = NUM_ESSENTIALS; + for (const auto &tot_name : solution_order) { + auto it = solution->Get_totals().find(tot_name); + if (it == solution->Get_totals().end()) { + data[i++] = 0.0; + continue; + } + data[i++] = it->second; + } +} + +void SolutionWrapper::set(const std::span &data) { + std::size_t i = NUM_ESSENTIALS; + cxxNameDouble new_totals; + for (const auto &tot_name : solution_order) { + const double value = data[i++]; + + if (value < 1E-25) { + continue; + } + new_totals[tot_name] = value; + } + + this->solution->Update(data[0], data[1], data[2], new_totals); +} + +std::vector +SolutionWrapper::names(cxxSolution *solution, + std::vector &solution_order) { + std::vector names; + + names.push_back("H"); + names.push_back("O"); + names.push_back("Charge"); + + for (const auto &[tot_name, value] : solution->Get_totals()) { + if (tot_name == "H(0)" || tot_name == "O(0)") { + continue; + } + + solution_order.push_back(tot_name); + } + + names.insert(names.end(), solution_order.begin(), solution_order.end()); + return names; +} \ No newline at end of file diff --git a/poet/src/SurfaceChargeWrapper.cpp b/poet/src/SurfaceChargeWrapper.cpp new file mode 100644 index 00000000..1eccfe7e --- /dev/null +++ b/poet/src/SurfaceChargeWrapper.cpp @@ -0,0 +1,72 @@ +#include "SurfaceWrapper.hpp" + +SurfaceWrapper::SurfaceChargeWrapper::SurfaceChargeWrapper( + cxxSurfaceCharge &charge, const std::set &solution_primaries) + : surface_charge(charge), primaries(solution_primaries) { + num_elements = NUM_NOT_TOTALS + solution_primaries.size(); +} + +void SurfaceWrapper::SurfaceChargeWrapper::get( + std::span &surface) const { + surface[0] = this->surface_charge.Get_specific_area(); + surface[1] = this->surface_charge.Get_grams(); + surface[2] = this->surface_charge.Get_charge_balance(); + surface[3] = this->surface_charge.Get_mass_water(); + surface[4] = this->surface_charge.Get_la_psi(); + + const auto dl_map = this->surface_charge.Get_diffuse_layer_totals(); + + for (auto it = this->primaries.begin(); it != this->primaries.end(); it++) { + const std::size_t index = + NUM_NOT_TOTALS + std::distance(this->primaries.begin(), it); + + const auto dl_it = dl_map.find(*it); + + if (dl_it == dl_map.end()) { + surface[index] = 0; + continue; + } + surface[index] = dl_it->second; + } +} + +void SurfaceWrapper::SurfaceChargeWrapper::set( + const std::span &surface) { + this->surface_charge.Set_specific_area(surface[0]); + this->surface_charge.Set_grams(surface[1]); + this->surface_charge.Set_charge_balance(surface[2]); + this->surface_charge.Set_mass_water(surface[3]); + this->surface_charge.Set_la_psi(surface[4]); + + auto &dl_map = this->surface_charge.Get_diffuse_layer_totals(); + dl_map.clear(); + + for (auto it = this->primaries.begin(); it != this->primaries.end(); it++) { + const std::size_t index = + NUM_NOT_TOTALS + std::distance(this->primaries.begin(), it); + + if (surface[index] == 0) { + continue; + } + + dl_map[*it] = surface[index]; + } +} + +std::vector SurfaceWrapper::SurfaceChargeWrapper::names( + cxxSurfaceCharge *charge, const std::set &primaries) { + std::vector names; + + const std::string &charge_name = charge->Get_name(); + names.push_back(charge_name + "_area"); + names.push_back(charge_name + "_grams"); + names.push_back(charge_name + "_cb"); + names.push_back(charge_name + "_mw"); + names.push_back(charge_name + "_la"); + + for (const auto &name : primaries) { + names.push_back(charge_name + "_tot_" + name); + } + + return names; +} \ No newline at end of file diff --git a/poet/src/SurfaceCompWrapper.cpp b/poet/src/SurfaceCompWrapper.cpp new file mode 100644 index 00000000..2bc45d1d --- /dev/null +++ b/poet/src/SurfaceCompWrapper.cpp @@ -0,0 +1,52 @@ +#include "SurfaceWrapper.hpp" + +SurfaceWrapper::SurfaceCompWrapper::SurfaceCompWrapper(cxxSurfaceComp &comp) + : surface_comp(comp) { + num_elements = NUM_NOT_TOTALS + comp.Get_totals().size(); + + for (const auto &[name, value] : comp.Get_totals()) { + total_names.push_back(name); + } +} + +void SurfaceWrapper::SurfaceCompWrapper::get(std::span &surface) const { + surface[0] = this->surface_comp.Get_moles(); + surface[1] = this->surface_comp.Get_la(); + surface[2] = this->surface_comp.Get_charge_balance(); + + const auto &totals = this->surface_comp.Get_totals(); + for (std::size_t i = 0; i < this->total_names.size(); i++) { + surface[NUM_NOT_TOTALS + i] = totals.at(this->total_names[i]); + } +} + +void SurfaceWrapper::SurfaceCompWrapper::set(const std::span &surface) { + this->surface_comp.Set_moles(surface[0]); + this->surface_comp.Set_la(surface[1]); + this->surface_comp.Set_charge_balance(surface[2]); + + auto &totals = this->surface_comp.Get_totals(); + totals.clear(); + + for (std::size_t i = 0; i < this->total_names.size(); i++) { + const std::size_t index = NUM_NOT_TOTALS + i; + + totals[this->total_names[i]] = surface[index]; + } +} + +std::vector +SurfaceWrapper::SurfaceCompWrapper::names(cxxSurfaceComp &comp) { + std::vector names; + + const std::string &phase_name = comp.Get_formula(); + names.push_back(phase_name + "_moles"); + names.push_back(phase_name + "_la"); + names.push_back(phase_name + "_cb"); + + for (const auto &tot : comp.Get_totals()) { + names.push_back(phase_name + "_" + tot.first); + } + + return names; +} \ No newline at end of file diff --git a/poet/src/SurfaceWrapper.cpp b/poet/src/SurfaceWrapper.cpp new file mode 100644 index 00000000..6fe396a9 --- /dev/null +++ b/poet/src/SurfaceWrapper.cpp @@ -0,0 +1,98 @@ +#include "SurfaceWrapper.hpp" +#include "SurfaceComp.h" + +SurfaceWrapper::SurfaceWrapper(cxxSurface *surf, + const std::set &solution_primaries, + const std::vector &comp_formulas, + const std::vector &charge_names) + : surface(surf) { + + for (const auto &comp_name : comp_formulas) { + auto it = std::find_if(surface->Get_surface_comps().begin(), + surface->Get_surface_comps().end(), + [&](const cxxSurfaceComp &comp) { + return comp.Get_formula() == comp_name; + }); + + if (it == surface->Get_surface_comps().end()) { + throw std::runtime_error( + "Surface component not found in Phreeqc variables"); + } + + surface_comps.push_back(std::make_unique(*it)); + + num_elements += surface_comps.back()->size(); + } + + for (const auto &charge_name : charge_names) { + auto it = std::find_if(surface->Get_surface_charges().begin(), + surface->Get_surface_charges().end(), + [&](const cxxSurfaceCharge &charge) { + return charge.Get_name() == charge_name; + }); + + if (it == surface->Get_surface_charges().end()) { + throw std::runtime_error("Surface charge not found in Phreeqc variables"); + } + + surface_charges.push_back( + std::make_unique(*it, solution_primaries)); + + num_elements += surface_charges.back()->size(); + } +} + +void SurfaceWrapper::get(std::span &surface) const { + std::size_t offset = 0; + + for (const auto &comp : surface_comps) { + std::span comp_span = surface.subspan(offset, comp->size()); + comp->get(comp_span); + offset += comp->size(); + } + + for (const auto &charge : surface_charges) { + std::span charge_span = surface.subspan(offset, charge->size()); + charge->get(charge_span); + offset += charge->size(); + } +} + +void SurfaceWrapper::set(const std::span &surface) { + std::size_t offset = 0; + + for (const auto &comp : surface_comps) { + std::span comp_span = surface.subspan(offset, comp->size()); + comp->set(comp_span); + offset += comp->size(); + } + + for (const auto &charge : surface_charges) { + std::span charge_span = surface.subspan(offset, charge->size()); + charge->set(charge_span); + offset += charge->size(); + } +} + +std::vector +SurfaceWrapper::names(cxxSurface *surface, + const std::set &solution_primaries, + std::vector &comp_formulas, + std::vector &charge_names) { + std::vector names; + + for (auto &comp : surface->Get_surface_comps()) { + comp_formulas.push_back(comp.Get_formula()); + const auto charge_names = SurfaceCompWrapper::names(comp); + names.insert(names.end(), charge_names.begin(), charge_names.end()); + } + + for (auto &charge : surface->Get_surface_charges()) { + charge_names.push_back(charge.Get_name()); + const auto comp_names = + SurfaceChargeWrapper::names(&charge, solution_primaries); + names.insert(names.end(), comp_names.begin(), comp_names.end()); + } + + return names; +} \ No newline at end of file diff --git a/src/phreeqcpp/Exchange.cxx b/src/phreeqcpp/Exchange.cxx index eae13837..7774ddf3 100644 --- a/src/phreeqcpp/Exchange.cxx +++ b/src/phreeqcpp/Exchange.cxx @@ -2,17 +2,18 @@ // ////////////////////////////////////////////////////////////////////// #ifdef _DEBUG -#pragma warning(disable : 4786) // disable truncation warning (Only used by debugger) +#pragma warning( \ + disable : 4786) // disable truncation warning (Only used by debugger) #endif -#include // std::cout std::cerr -#include // assert -#include // std::sort +#include // std::sort +#include // assert +#include // std::cout std::cerr #include "Utils.h" // define first +#include "Exchange.h" #include "Phreeqc.h" #include "cxxMix.h" -#include "Exchange.h" #include "phqalloc.h" #if defined(PHREEQCI_GUI) @@ -28,493 +29,367 @@ static char THIS_FILE[] = __FILE__; ////////////////////////////////////////////////////////////////////// cxxExchange::cxxExchange(PHRQ_io *io) - // - // default constructor for cxxExchange - // -: cxxNumKeyword(io) -{ - new_def = false; - solution_equilibria = false; - n_solution = -999; - pitzer_exchange_gammas = true; + // + // default constructor for cxxExchange + // + : cxxNumKeyword(io) { + new_def = false; + solution_equilibria = false; + n_solution = -999; + pitzer_exchange_gammas = true; } -cxxExchange::cxxExchange(const std::map < int, cxxExchange > &entities, - cxxMix & mix, int l_n_user, PHRQ_io *io): -cxxNumKeyword(io) -{ - this->n_user = this->n_user_end = l_n_user; - this->pitzer_exchange_gammas = true; - this->new_def = false; - this->n_solution = -999; - this->solution_equilibria = false; -// -// Mix exchangers -// - const std::map < int, LDBLE >&mixcomps = mix.Get_mixComps(); - std::map < int, LDBLE >::const_iterator it; - for (it = mixcomps.begin(); it != mixcomps.end(); it++) - { - if (entities.find(it->first) != entities.end()) - { - const cxxExchange *entity_ptr = - &(entities.find(it->first)->second); - this->add(*entity_ptr, it->second); - this->pitzer_exchange_gammas = entity_ptr->pitzer_exchange_gammas; - } - } -} - - -cxxExchange::~cxxExchange() -{ -} -bool -cxxExchange::Get_related_phases() const -{ - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - if (this->exchange_comps[i].Get_phase_name().size() > 0) - { - return (true); - } - } - return (false); -} - -bool -cxxExchange::Get_related_rate() const -{ - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - if (this->exchange_comps[i].Get_rate_name().size() > 0) - { - return (true); - } - } - return (false); -} - -void -cxxExchange::dump_xml(std::ostream & s_oss, unsigned int indent) const -{ - unsigned int i; - s_oss.precision(DBL_DIG - 1); - std::string indent0(""), indent1(""), indent2(""); - for (i = 0; i < indent; ++i) - indent0.append(Utilities::INDENT); - for (i = 0; i < indent + 1; ++i) - indent1.append(Utilities::INDENT); - for (i = 0; i < indent + 2; ++i) - indent2.append(Utilities::INDENT); - - // Exchange element and attributes - s_oss << indent0; - s_oss << " - pitzer_exchange_gammas << "\"" << "\n"; - - // components - s_oss << indent1; - s_oss << "exchange_comps.size(); j++) - { - this->exchange_comps[j].dump_xml(s_oss, indent + 2); - } - - return; -} -void -cxxExchange::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) const -{ - unsigned int i; - s_oss.precision(DBL_DIG - 1); - std::string indent0(""), indent1(""), indent2(""); - for (i = 0; i < indent; ++i) - indent0.append(Utilities::INDENT); - for (i = 0; i < indent + 1; ++i) - indent1.append(Utilities::INDENT); - for (i = 0; i < indent + 2; ++i) - indent2.append(Utilities::INDENT); - - // Exchange element and attributes - s_oss << indent0; - int n_user_local = (n_out != NULL) ? *n_out : this->n_user; - s_oss << "EXCHANGE_RAW " << n_user_local << " " << this->description << "\n"; - - s_oss << indent1 << "# EXCHANGE_MODIFY candidate identifiers #\n"; - s_oss << indent1; - s_oss << "-exchange_gammas " << (this->pitzer_exchange_gammas ? 1 : 0) << "\n"; - // exchComps - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - s_oss << indent1; - s_oss << "-component " << this->exchange_comps[i].Get_formula() << "\n"; - this->exchange_comps[i].dump_raw(s_oss, indent + 2); - } - - s_oss << indent1 << "# EXCHANGE_MODIFY candidates with new_def=true #\n"; - s_oss << indent1; - s_oss << "-new_def " << 0 << "\n"; - s_oss << indent1; - s_oss << "-solution_equilibria " << 0 << "\n"; - s_oss << indent1; - s_oss << "-n_solution " << this->n_solution << "\n"; - - s_oss << indent1 << "# Exchange workspace variables #\n"; - s_oss << indent1; - s_oss << "-totals" << "\n"; - this->totals.dump_raw(s_oss, indent + 1); - - return; -} - -void cxxExchange::dump_essential_names( - std::vector &e_names) const { - e_names.clear(); - e_names.reserve(this->exchange_comps.size() * - this->exchange_comps.at(0).Get_totals().size()); - - for (const auto &comp : this->exchange_comps) { - const std::string formular = comp.Get_formula(); - e_names.push_back(formular); - for (const auto &total : comp.Get_totals()) { - if (total.first == formular) { - continue; - } - e_names.push_back(total.first + formular); +cxxExchange::cxxExchange(const std::map &entities, + cxxMix &mix, int l_n_user, PHRQ_io *io) + : cxxNumKeyword(io) { + this->n_user = this->n_user_end = l_n_user; + this->pitzer_exchange_gammas = true; + this->new_def = false; + this->n_solution = -999; + this->solution_equilibria = false; + // + // Mix exchangers + // + const std::map &mixcomps = mix.Get_mixComps(); + std::map::const_iterator it; + for (it = mixcomps.begin(); it != mixcomps.end(); it++) { + if (entities.find(it->first) != entities.end()) { + const cxxExchange *entity_ptr = &(entities.find(it->first)->second); + this->add(*entity_ptr, it->second); + this->pitzer_exchange_gammas = entity_ptr->pitzer_exchange_gammas; } } } -void cxxExchange::get_essential_values(std::vector &e_values) const { - e_values.clear(); - e_values.reserve(this->exchange_comps.size() * - this->exchange_comps.at(0).Get_totals().size()); - - for (const auto &comp : this->exchange_comps) { - const std::string formular = comp.Get_formula(); - e_values.push_back(comp.Get_totals().find(formular)->second); - for (const auto &total : comp.Get_totals()) { - if (total.first == formular) { - continue; - } - e_values.push_back(total.second); +cxxExchange::~cxxExchange() {} +bool cxxExchange::Get_related_phases() const { + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + if (this->exchange_comps[i].Get_phase_name().size() > 0) { + return (true); } } + return (false); +} + +bool cxxExchange::Get_related_rate() const { + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + if (this->exchange_comps[i].Get_rate_name().size() > 0) { + return (true); + } + } + return (false); +} + +void cxxExchange::dump_xml(std::ostream &s_oss, unsigned int indent) const { + unsigned int i; + s_oss.precision(DBL_DIG - 1); + std::string indent0(""), indent1(""), indent2(""); + for (i = 0; i < indent; ++i) + indent0.append(Utilities::INDENT); + for (i = 0; i < indent + 1; ++i) + indent1.append(Utilities::INDENT); + for (i = 0; i < indent + 2; ++i) + indent2.append(Utilities::INDENT); + + // Exchange element and attributes + s_oss << indent0; + s_oss << "pitzer_exchange_gammas << "\"" + << "\n"; + + // components + s_oss << indent1; + s_oss << "exchange_comps.size(); j++) { + this->exchange_comps[j].dump_xml(s_oss, indent + 2); + } return; } +void cxxExchange::dump_raw(std::ostream &s_oss, unsigned int indent, + int *n_out) const { + unsigned int i; + s_oss.precision(DBL_DIG - 1); + std::string indent0(""), indent1(""), indent2(""); + for (i = 0; i < indent; ++i) + indent0.append(Utilities::INDENT); + for (i = 0; i < indent + 1; ++i) + indent1.append(Utilities::INDENT); + for (i = 0; i < indent + 2; ++i) + indent2.append(Utilities::INDENT); -void cxxExchange::set_essential_values(std::vector::iterator &it) { - for (auto &comp : this->exchange_comps) { - const std::string formular = comp.Get_formula(); - auto formular_total = comp.Get_totals(); - formular_total[formular] = *(it++); - for (auto &total : comp.Get_totals()) { - if (total.first == formular) { - continue; - } - total.second = *(it++); - } + // Exchange element and attributes + s_oss << indent0; + int n_user_local = (n_out != NULL) ? *n_out : this->n_user; + s_oss << "EXCHANGE_RAW " << n_user_local << " " + << this->description << "\n"; + + s_oss << indent1 << "# EXCHANGE_MODIFY candidate identifiers #\n"; + s_oss << indent1; + s_oss << "-exchange_gammas " + << (this->pitzer_exchange_gammas ? 1 : 0) << "\n"; + // exchComps + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + s_oss << indent1; + s_oss << "-component " + << this->exchange_comps[i].Get_formula() << "\n"; + this->exchange_comps[i].dump_raw(s_oss, indent + 2); } + s_oss << indent1 << "# EXCHANGE_MODIFY candidates with new_def=true #\n"; + s_oss << indent1; + s_oss << "-new_def " << 0 << "\n"; + s_oss << indent1; + s_oss << "-solution_equilibria " << 0 << "\n"; + s_oss << indent1; + s_oss << "-n_solution " << this->n_solution << "\n"; + + s_oss << indent1 << "# Exchange workspace variables #\n"; + s_oss << indent1; + s_oss << "-totals" + << "\n"; + this->totals.dump_raw(s_oss, indent + 1); + return; } -void -cxxExchange::read_raw(CParser & parser, bool check) -{ - std::istream::pos_type ptr; - std::istream::pos_type next_char; - std::string token; - bool useLastLine(false); +void cxxExchange::read_raw(CParser &parser, bool check) { + std::istream::pos_type ptr; + std::istream::pos_type next_char; + std::string token; + bool useLastLine(false); - // Read exchange number and description - this->read_number_description(parser); - this->Set_new_def(false); - - bool pitzer_exchange_gammas_defined(false); + // Read exchange number and description + this->read_number_description(parser); + this->Set_new_def(false); - for (;;) - { - int opt; - if (useLastLine == false) - { - opt = parser.get_option(vopts, next_char); - } - else - { - opt = parser.getOptionFromLastLine(vopts, next_char, true); - } - useLastLine = false; - switch (opt) - { - case CParser::OPT_EOF: - break; - case CParser::OPT_KEYWORD: - break; - case CParser::OPT_DEFAULT: - case CParser::OPT_ERROR: - opt = CParser::OPT_EOF; - parser.error_msg("Unknown input in EXCH_COMP_RAW keyword.", - PHRQ_io::OT_CONTINUE); - parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE); - break; + bool pitzer_exchange_gammas_defined(false); - case 0: // pitzer_exchange_gammas - case 2: // exchange_gammas - if (!(parser.get_iss() >> this->pitzer_exchange_gammas)) - { - this->pitzer_exchange_gammas = false; - parser.incr_input_error(); - parser. - error_msg - ("Expected boolean value for pitzer_exchange_gammas.", - PHRQ_io::OT_CONTINUE); - } - pitzer_exchange_gammas_defined = true; - break; - case 1: // component - { - std::string str; - if (!(parser.get_iss() >> str)) - { - parser.incr_input_error(); - parser.error_msg("Expected string value for component name.", - PHRQ_io::OT_CONTINUE); - } - else - { - cxxExchComp temp_comp(this->io); - temp_comp.Set_formula(str.c_str()); - cxxExchComp *comp_ptr = this->Find_comp(str); - if (comp_ptr) - { - temp_comp = *comp_ptr; - } - temp_comp.read_raw(parser, check); - if (comp_ptr) - { - for (size_t j = 0; j < this->exchange_comps.size(); j++) - { - if (Utilities::strcmp_nocase(this->exchange_comps[j].Get_formula().c_str(), str.c_str()) == 0) - { - this->exchange_comps[j] = temp_comp; - } - } - } - else - { - this->exchange_comps.push_back(temp_comp); - } - useLastLine = true; - } - } - break; - case 3: // new_def - if (!(parser.get_iss() >> this->new_def)) - { - this->new_def = false; - parser.incr_input_error(); - parser. - error_msg - ("Expected boolean value for new_def.", - PHRQ_io::OT_CONTINUE); - } - break; - case 4: // solution_equilibria - if (!(parser.get_iss() >> this->solution_equilibria)) - { - this->solution_equilibria = false; - parser.incr_input_error(); - parser. - error_msg - ("Expected boolean value for solution_equilibria.", - PHRQ_io::OT_CONTINUE); - } - break; - case 5: // n_solution - if (!(parser.get_iss() >> this->n_solution)) - { - this->n_solution = -999; - parser.incr_input_error(); - parser. - error_msg - ("Expected integer value for n_solution.", - PHRQ_io::OT_CONTINUE); - } - break; - case 6: // totals - if (this->totals.read_raw(parser, next_char) != - CParser::PARSER_OK) - { - parser.incr_input_error(); - parser. - error_msg - ("Expected element name and molality for Exchange totals.", - PHRQ_io::OT_CONTINUE); - } - break; - } - if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) - break; - } - // members that must be defined - if (check) - { - if (pitzer_exchange_gammas_defined == false) - { - parser.incr_input_error(); - parser. - error_msg - ("Pitzer_exchange_gammsa not defined for EXCHANGE_RAW input.", - PHRQ_io::OT_CONTINUE); - } - } - this->Sort_comps(); + for (;;) { + int opt; + if (useLastLine == false) { + opt = parser.get_option(vopts, next_char); + } else { + opt = parser.getOptionFromLastLine(vopts, next_char, true); + } + useLastLine = false; + switch (opt) { + case CParser::OPT_EOF: + break; + case CParser::OPT_KEYWORD: + break; + case CParser::OPT_DEFAULT: + case CParser::OPT_ERROR: + opt = CParser::OPT_EOF; + parser.error_msg("Unknown input in EXCH_COMP_RAW keyword.", + PHRQ_io::OT_CONTINUE); + parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE); + break; + + case 0: // pitzer_exchange_gammas + case 2: // exchange_gammas + if (!(parser.get_iss() >> this->pitzer_exchange_gammas)) { + this->pitzer_exchange_gammas = false; + parser.incr_input_error(); + parser.error_msg("Expected boolean value for pitzer_exchange_gammas.", + PHRQ_io::OT_CONTINUE); + } + pitzer_exchange_gammas_defined = true; + break; + case 1: // component + { + std::string str; + if (!(parser.get_iss() >> str)) { + parser.incr_input_error(); + parser.error_msg("Expected string value for component name.", + PHRQ_io::OT_CONTINUE); + } else { + cxxExchComp temp_comp(this->io); + temp_comp.Set_formula(str.c_str()); + cxxExchComp *comp_ptr = this->Find_comp(str); + if (comp_ptr) { + temp_comp = *comp_ptr; + } + temp_comp.read_raw(parser, check); + if (comp_ptr) { + for (size_t j = 0; j < this->exchange_comps.size(); j++) { + if (Utilities::strcmp_nocase( + this->exchange_comps[j].Get_formula().c_str(), + str.c_str()) == 0) { + this->exchange_comps[j] = temp_comp; + } + } + } else { + this->exchange_comps.push_back(temp_comp); + } + useLastLine = true; + } + } break; + case 3: // new_def + if (!(parser.get_iss() >> this->new_def)) { + this->new_def = false; + parser.incr_input_error(); + parser.error_msg("Expected boolean value for new_def.", + PHRQ_io::OT_CONTINUE); + } + break; + case 4: // solution_equilibria + if (!(parser.get_iss() >> this->solution_equilibria)) { + this->solution_equilibria = false; + parser.incr_input_error(); + parser.error_msg("Expected boolean value for solution_equilibria.", + PHRQ_io::OT_CONTINUE); + } + break; + case 5: // n_solution + if (!(parser.get_iss() >> this->n_solution)) { + this->n_solution = -999; + parser.incr_input_error(); + parser.error_msg("Expected integer value for n_solution.", + PHRQ_io::OT_CONTINUE); + } + break; + case 6: // totals + if (this->totals.read_raw(parser, next_char) != CParser::PARSER_OK) { + parser.incr_input_error(); + parser.error_msg( + "Expected element name and molality for Exchange totals.", + PHRQ_io::OT_CONTINUE); + } + break; + } + if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) + break; + } + // members that must be defined + if (check) { + if (pitzer_exchange_gammas_defined == false) { + parser.incr_input_error(); + parser.error_msg( + "Pitzer_exchange_gammsa not defined for EXCHANGE_RAW input.", + PHRQ_io::OT_CONTINUE); + } + } + this->Sort_comps(); } -void -cxxExchange::add(const cxxExchange & addee, LDBLE extensive) - // - // Add existing exchange to "this" exchange - // +void cxxExchange::add(const cxxExchange &addee, LDBLE extensive) +// +// Add existing exchange to "this" exchange +// { - // exchComps - if (extensive == 0.0) - return; - for (size_t i = 0; i < addee.exchange_comps.size(); i++) - { - size_t j; - for (j = 0; j < this->Get_exchange_comps().size(); j++) - if (addee.exchange_comps[i].Get_formula() == this->exchange_comps[j].Get_formula()) - { - this->exchange_comps[j].add(addee.exchange_comps[i], extensive); - break; - } - if (j == this->exchange_comps.size()) - { - cxxExchComp exc = addee.exchange_comps[i]; - exc.multiply(extensive); - this->exchange_comps.push_back(exc); - } - } - this->pitzer_exchange_gammas = addee.pitzer_exchange_gammas; + // exchComps + if (extensive == 0.0) + return; + for (size_t i = 0; i < addee.exchange_comps.size(); i++) { + size_t j; + for (j = 0; j < this->Get_exchange_comps().size(); j++) + if (addee.exchange_comps[i].Get_formula() == + this->exchange_comps[j].Get_formula()) { + this->exchange_comps[j].add(addee.exchange_comps[i], extensive); + break; + } + if (j == this->exchange_comps.size()) { + cxxExchComp exc = addee.exchange_comps[i]; + exc.multiply(extensive); + this->exchange_comps.push_back(exc); + } + } + this->pitzer_exchange_gammas = addee.pitzer_exchange_gammas; } -void -cxxExchange::totalize() -{ - this->totals.clear(); - // component structures - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - this->totals.add_extensive(this->exchange_comps[i].Get_totals(), 1.0); - this->totals.add("Charge", this->exchange_comps[i].Get_charge_balance()); - } - return; +void cxxExchange::totalize() { + this->totals.clear(); + // component structures + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + this->totals.add_extensive(this->exchange_comps[i].Get_totals(), 1.0); + this->totals.add("Charge", this->exchange_comps[i].Get_charge_balance()); + } + return; } -bool -cxxExchange::Get_pitzer_exchange_gammas() const -{ - return this->pitzer_exchange_gammas; +bool cxxExchange::Get_pitzer_exchange_gammas() const { + return this->pitzer_exchange_gammas; } -void -cxxExchange::Set_pitzer_exchange_gammas(bool b) -{ - this->pitzer_exchange_gammas = b; +void cxxExchange::Set_pitzer_exchange_gammas(bool b) { + this->pitzer_exchange_gammas = b; } -const cxxNameDouble & -cxxExchange::Get_totals() const -{ - return totals; +const cxxNameDouble &cxxExchange::Get_totals() const { return totals; } +cxxExchComp *cxxExchange::Find_comp(std::string s) { + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + cxxNameDouble nd(this->exchange_comps[i].Get_totals()); + cxxNameDouble::iterator nd_it; + for (nd_it = nd.begin(); nd_it != nd.end(); nd_it++) { + if (nd_it->first == s) { + return (&(this->exchange_comps[i])); + } + } + } + return NULL; } -cxxExchComp *cxxExchange::Find_comp(std::string s) -{ - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - cxxNameDouble nd(this->exchange_comps[i].Get_totals()); - cxxNameDouble::iterator nd_it; - for (nd_it = nd.begin(); nd_it != nd.end(); nd_it++) - { - if(nd_it->first == s) - { - return (&(this->exchange_comps[i])); - } - } - } - return NULL; -} -void cxxExchange:: -Sort_comps(void) -{ - // sort comps - { - std::map comp_map; - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - comp_map[this->exchange_comps[i].Get_formula()] = this->exchange_comps[i]; - } - this->exchange_comps.clear(); - std::map::iterator it; - for (it = comp_map.begin(); it != comp_map.end(); it++) - { - this->exchange_comps.push_back(it->second); - } - } +void cxxExchange::Sort_comps(void) { + // sort comps + { + std::map comp_map; + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + comp_map[this->exchange_comps[i].Get_formula()] = this->exchange_comps[i]; + } + this->exchange_comps.clear(); + std::map::iterator it; + for (it = comp_map.begin(); it != comp_map.end(); it++) { + this->exchange_comps.push_back(it->second); + } + } } /* ---------------------------------------------------------------------- */ -void -cxxExchange::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles) +void cxxExchange::Serialize(Dictionary &dictionary, std::vector &ints, + std::vector &doubles) /* ---------------------------------------------------------------------- */ { - ints.push_back(this->n_user); - ints.push_back((int) this->exchange_comps.size()); - for (size_t i = 0; i < this->exchange_comps.size(); i++) - { - exchange_comps[i].Serialize(dictionary, ints, doubles); - } - ints.push_back(this->pitzer_exchange_gammas ? 1 : 0); - ints.push_back(this->new_def ? 1 : 0); - ints.push_back(this->solution_equilibria ? 1 : 0); - ints.push_back(this->n_solution); - this->totals.Serialize(dictionary, ints, doubles); - + ints.push_back(this->n_user); + ints.push_back((int)this->exchange_comps.size()); + for (size_t i = 0; i < this->exchange_comps.size(); i++) { + exchange_comps[i].Serialize(dictionary, ints, doubles); + } + ints.push_back(this->pitzer_exchange_gammas ? 1 : 0); + ints.push_back(this->new_def ? 1 : 0); + ints.push_back(this->solution_equilibria ? 1 : 0); + ints.push_back(this->n_solution); + this->totals.Serialize(dictionary, ints, doubles); } /* ---------------------------------------------------------------------- */ -void -cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles, int &ii, int &dd) +void cxxExchange::Deserialize(Dictionary &dictionary, std::vector &ints, + std::vector &doubles, int &ii, int &dd) /* ---------------------------------------------------------------------- */ { - this->n_user = ints[ii++]; - this->n_user_end = this->n_user; - this->description = " "; - - int count = ints[ii++]; - this->exchange_comps.clear(); - for (int n = 0; n < count; n++) - { - cxxExchComp ec(this->io); - ec.Deserialize(dictionary, ints, doubles, ii, dd); - this->exchange_comps.push_back(ec); - } - this->pitzer_exchange_gammas = (ints[ii++] != 0); - this->new_def = (ints[ii++] != 0); - this->solution_equilibria = (ints[ii++] != 0); - this->n_solution = ints[ii++]; - this->totals.Deserialize(dictionary, ints, doubles, ii, dd); + this->n_user = ints[ii++]; + this->n_user_end = this->n_user; + this->description = " "; + int count = ints[ii++]; + this->exchange_comps.clear(); + for (int n = 0; n < count; n++) { + cxxExchComp ec(this->io); + ec.Deserialize(dictionary, ints, doubles, ii, dd); + this->exchange_comps.push_back(ec); + } + this->pitzer_exchange_gammas = (ints[ii++] != 0); + this->new_def = (ints[ii++] != 0); + this->solution_equilibria = (ints[ii++] != 0); + this->n_solution = ints[ii++]; + this->totals.Deserialize(dictionary, ints, doubles, ii, dd); } - -const std::vector< std::string >::value_type temp_vopts[] = { - std::vector< std::string >::value_type("pitzer_exchange_gammas"), // 0 - std::vector< std::string >::value_type("component"), // 1 - std::vector< std::string >::value_type("exchange_gammas"), // 2 - std::vector< std::string >::value_type("new_def"), // 3 - std::vector< std::string >::value_type("solution_equilibria"), // 4 - std::vector< std::string >::value_type("n_solution"), // 5 - std::vector< std::string >::value_type("totals") // 6 -}; -const std::vector< std::string > cxxExchange::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); +const std::vector::value_type temp_vopts[] = { + std::vector::value_type("pitzer_exchange_gammas"), // 0 + std::vector::value_type("component"), // 1 + std::vector::value_type("exchange_gammas"), // 2 + std::vector::value_type("new_def"), // 3 + std::vector::value_type("solution_equilibria"), // 4 + std::vector::value_type("n_solution"), // 5 + std::vector::value_type("totals") // 6 +}; +const std::vector + cxxExchange::vopts(temp_vopts, + temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/src/phreeqcpp/Exchange.h b/src/phreeqcpp/Exchange.h index 8fe9f96b..770db96c 100644 --- a/src/phreeqcpp/Exchange.h +++ b/src/phreeqcpp/Exchange.h @@ -20,11 +20,8 @@ public: cxxMix & mx, int n_user, PHRQ_io *io=NULL); ~cxxExchange(); - void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const; - void dump_essential_names(std::vector &e_names) const; - - void get_essential_values(std::vector &e_values) const; - void set_essential_values(std::vector::iterator &it); + void dump_raw(std::ostream &s_oss, unsigned int indent, + int *n_out = NULL) const; void read_raw(CParser & parser, bool check = true); diff --git a/src/phreeqcpp/PPassemblage.cxx b/src/phreeqcpp/PPassemblage.cxx index 0145c6a6..8b9b830e 100644 --- a/src/phreeqcpp/PPassemblage.cxx +++ b/src/phreeqcpp/PPassemblage.cxx @@ -131,40 +131,6 @@ cxxPPassemblage::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) this->assemblage_totals.dump_raw(s_oss, indent + 1); } -void cxxPPassemblage::dump_essential_names( - std::vector &e_names) const { - e_names.clear(); - e_names.reserve(this->pp_assemblage_comps.size() * 2); - - for (const auto &pp_comp : this->pp_assemblage_comps) { - e_names.push_back(pp_comp.first); - e_names.push_back(pp_comp.first + "_si"); - } - - return; -} - -void cxxPPassemblage::get_essential_values(std::vector &e_values) const { - e_values.clear(); - e_values.reserve(this->pp_assemblage_comps.size() * 2); - - for (const auto &pp_comp : this->pp_assemblage_comps) { - e_values.push_back(pp_comp.second.Get_moles()); - e_values.push_back(pp_comp.second.Get_si()); - } - - return; -} - -void cxxPPassemblage::set_essential_values(std::vector::iterator &it) { - for (auto &pp_comp : this->pp_assemblage_comps) { - pp_comp.second.Set_moles(*(it++)); - pp_comp.second.Set_si(*(it++)); - } - - return; -} - void cxxPPassemblage::read_raw(CParser & parser, bool check) { diff --git a/src/phreeqcpp/PPassemblage.h b/src/phreeqcpp/PPassemblage.h index e19cb1ab..b9fcb214 100644 --- a/src/phreeqcpp/PPassemblage.h +++ b/src/phreeqcpp/PPassemblage.h @@ -23,10 +23,6 @@ class cxxPPassemblage:public cxxNumKeyword void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const; - void dump_essential_names(std::vector &e_names) const; - void get_essential_values(std::vector &e_values) const; - void set_essential_values(std::vector::iterator &it); - void read_raw(CParser &parser, bool check = true); const cxxNameDouble & Get_assemblage_totals() const diff --git a/src/phreeqcpp/Phreeqc.h b/src/phreeqcpp/Phreeqc.h index 16b196c9..79f16658 100644 --- a/src/phreeqcpp/Phreeqc.h +++ b/src/phreeqcpp/Phreeqc.h @@ -424,7 +424,10 @@ public: int initial_exchangers(int print); int initial_gas_phases(int print); int initial_solutions(int print); - int step_save_exch(int n_user); + + int initial_solutions_poet(int sol_id); + + int step_save_exch(int n_user); int step_save_surf(int n_user); int initial_surfaces(int print); int reactions(void); @@ -1830,8 +1833,10 @@ protected: friend class IPhreeqcMMS; friend class IPhreeqcPhast; friend class PhreeqcRM; + friend class PhreeqcInit; + friend class PhreeqcEngine; - std::vector keycount; // used to mark keywords that have been read + std::vector keycount; // used to mark keywords that have been read public: static const class const_iso iso_defaults[]; diff --git a/src/phreeqcpp/Solution.cxx b/src/phreeqcpp/Solution.cxx index 66c7f29c..01bc6053 100644 --- a/src/phreeqcpp/Solution.cxx +++ b/src/phreeqcpp/Solution.cxx @@ -367,70 +367,6 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con return; } -void cxxSolution::dump_essential_names( - std::vector &e_names) const { - e_names.clear(); - e_names.reserve(this->totals.size() + 3); - - // e_names.push_back("H20"); - e_names.push_back("H"); - e_names.push_back("O"); - e_names.push_back("Charge"); - - for (const auto &total : this->totals) { - if ((total.first == "O(0)") || (total.first == "H(0)")) { - continue; - } - e_names.push_back(total.first); - } - - return; -} - -void cxxSolution::get_essential_values( - std::vector &e_values, const std::vector &order) const { - e_values.clear(); - e_values.reserve(this->totals.size() + 4); - - // e_values.push_back(mass_water); - e_values.push_back(total_h); - e_values.push_back(total_o); - e_values.push_back(cb); - - for (int i = 3; i < order.size(); i++) { - const auto tot = totals.find(order[i]); - if (tot == totals.end()) { - e_values.push_back(0); - continue; - } - e_values.push_back((tot->second)); - } - - return; -} - -void cxxSolution::set_essential_values(std::vector::iterator &it, - const std::vector &order) { - // skip first field of moles H2O - // it++; - - double t_h = *(it++); - double t_o = *(it++); - double c_b = *(it++); - - cxxNameDouble nd; - for (int i = 3; i < order.size(); i++) { - double curr_val = *(it++); - - if (curr_val < 0) { - curr_val = 0; - } - nd.add(order[i].c_str(), curr_val); - } - - this->Update(t_h, t_o, c_b, nd); -} - #ifdef USE_REVISED_READ_RAW void cxxSolution::read_raw(CParser & parser, bool check) @@ -1258,17 +1194,20 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble this->cb = charge; this->mass_water = o_tot / 55.5; - // Don`t bother to update activities? - this->Update(const_nd); - //this->totals = const_nd; - cxxNameDouble::iterator it; - for (it = this->totals.begin(); it != this->totals.end(); it++) - { - if (it->second < 1e-25) - { - it->second = 0.0; - } - } + // Don`t bother to update activities? + // this->Update(const_nd); + + // UP: instead of updating, we just set the new totals ... + this->totals = const_nd; + + // UP: ... and doesn't need another check for zero values + + // cxxNameDouble::iterator it; + // for (it = this->totals.begin(); it != this->totals.end(); it++) { + // if (it->second < 1e-25) { + // it->second = 0.0; + // } + // } } void cxxSolution::Update_activities(const cxxNameDouble &original_tot) diff --git a/src/phreeqcpp/Solution.h b/src/phreeqcpp/Solution.h index c821f64f..7bb1f740 100644 --- a/src/phreeqcpp/Solution.h +++ b/src/phreeqcpp/Solution.h @@ -105,13 +105,6 @@ class cxxSolution:public cxxNumKeyword void dump_xml(std::ostream & os, unsigned int indent = 0) const; void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const; - void dump_essential_names(std::vector &e_names) const; - - void get_essential_values(std::vector &e_values, - const std::vector &order) const; - void set_essential_values(std::vector::iterator &it, - const std::vector &order); - LDBLE Get_master_activity(char *string) const; void Set_master_activity(char *string, LDBLE value); LDBLE Get_total(const char *string) const; diff --git a/src/phreeqcpp/Surface.cxx b/src/phreeqcpp/Surface.cxx index 94a310bf..ed55bb2d 100644 --- a/src/phreeqcpp/Surface.cxx +++ b/src/phreeqcpp/Surface.cxx @@ -178,105 +178,6 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons return; } -void cxxSurface::dump_essential_names(std::vector &e_names) { - // this is just some reserve size, might be changed in the future - const int reserve_size = - surface_comps.size() * 4 + surface_charges.size() * 100; - e_names.clear(); - e_names.reserve(reserve_size); - - for (auto &comp : this->surface_comps) { - const std::string phase_name = comp.Get_formula(); - e_names.push_back(phase_name + "_moles"); - e_names.push_back(phase_name + "_la"); - e_names.push_back(phase_name + "_cb"); - - for (const auto &tot : comp.Get_totals()) { - e_names.push_back(phase_name + "_" + tot.first); - } - } - - for (auto &charge : this->surface_charges) { - const std::string component = charge.Get_name() + "_charge"; - e_names.push_back(component + "_area"); - e_names.push_back(component + "_grams"); - e_names.push_back(component + "_cb"); - e_names.push_back(component + "_mw"); - e_names.push_back(component + "_la"); - - for (const auto &tot : charge.Get_diffuse_layer_totals()) { - e_names.push_back(component + "_tot_" + tot.first); - } - // for (const auto &dl_layer : charge.Get_dl_species_map()) { - // e_names.push_back(component + "_dls_" + - // std::to_string(dl_layer.first)); - // } - } -} - -void cxxSurface::get_essential_values(std::vector &e_values) { - // this is just some reserve size, might be changed in the future - const int reserve_size = - surface_comps.size() * 4 + surface_charges.size() * 100; - e_values.clear(); - e_values.reserve(reserve_size); - - for (auto &comp : this->surface_comps) { - e_values.push_back(comp.Get_moles()); - e_values.push_back(comp.Get_la()); - e_values.push_back(comp.Get_charge_balance()); - - for (const auto &tot : comp.Get_totals()) { - e_values.push_back(tot.second); - } - } - - for (auto &charge : this->surface_charges) { - e_values.push_back(charge.Get_specific_area()); - e_values.push_back(charge.Get_grams()); - e_values.push_back(charge.Get_charge_balance()); - e_values.push_back(charge.Get_mass_water()); - e_values.push_back(charge.Get_la_psi()); - - for (const auto &tot : charge.Get_diffuse_layer_totals()) { - e_values.push_back(tot.second); - } - // for (const auto &dl_layer : charge.Get_dl_species_map()) { - // e_values.push_back(dl_layer.second); - // } - } -} - -void cxxSurface::set_essential_values(std::vector::iterator &it) { - for (auto &comp : this->surface_comps) { - comp.Set_moles(*(it++)); - comp.Set_la(*(it++)); - comp.Set_charge_balance(*(it++)); - - for (auto &tot : comp.Get_totals()) { - tot.second = *(it++); - } - } - - for (auto &charge : this->surface_charges) { - charge.Set_specific_area(*(it++)); - charge.Set_grams(*(it++)); - charge.Set_charge_balance(*(it++)); - charge.Set_mass_water(*(it++)); - charge.Set_la_psi(*(it++)); - - for (auto &tot : charge.Get_diffuse_layer_totals()) { - tot.second = *(it++); - } - - charge.Get_dl_species_map().clear(); - - // for (auto &dl_layer : charge.Get_dl_species_map()) { - // dl_layer.second = *(it++); - // } - } -} - void cxxSurface::read_raw(CParser & parser, bool check) { diff --git a/src/phreeqcpp/Surface.h b/src/phreeqcpp/Surface.h index 0a61a3b5..04527265 100644 --- a/src/phreeqcpp/Surface.h +++ b/src/phreeqcpp/Surface.h @@ -1,3 +1,4 @@ +#include #if !defined(SURFACE_H_INCLUDED) #define SURFACE_H_INCLUDED @@ -31,11 +32,6 @@ public: //void dump_xml(std::ostream & os, unsigned int indent = 0) const; void dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out=NULL) const; - void dump_essential_names(std::vector &e_names); - - void get_essential_values(std::vector &e_values); - void set_essential_values(std::vector::iterator &it); - void read_raw(CParser & parser, bool check = true); bool Get_related_phases(void) const; bool Get_related_rate(void) const; diff --git a/src/phreeqcpp/cxxKinetics.cxx b/src/phreeqcpp/cxxKinetics.cxx index bdfb7d83..636111a4 100644 --- a/src/phreeqcpp/cxxKinetics.cxx +++ b/src/phreeqcpp/cxxKinetics.cxx @@ -159,45 +159,6 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out) co return; } -void cxxKinetics::dump_essential_names( - std::vector &e_names) const { - e_names.clear(); - e_names.reserve(this->kinetics_comps.size() * 2); - - for (const auto &comp : this->kinetics_comps) { - const std::string name = comp.Get_rate_name(); - e_names.push_back(name); - for (int i = 0; i < comp.Get_d_params().size(); i++) { - e_names.push_back(name + "_p" + std::to_string(i + 1)); - } - } -} - -void cxxKinetics::get_essential_values(std::vector &e_values) const { - e_values.clear(); - e_values.reserve(this->kinetics_comps.size() * 2); - - for (const auto &comp : this->kinetics_comps) { - e_values.push_back(comp.Get_m()); - for (const auto ¶m : comp.Get_d_params()) { - e_values.push_back(param); - } - } - - return; -} - -void cxxKinetics::set_essential_values(std::vector::iterator &it) { - for (auto &comp : this->kinetics_comps) { - comp.Set_m(*(it++)); - for (auto ¶m : comp.Get_d_params()) { - param = *(it++); - } - } - - return; -} - void cxxKinetics::read_raw(CParser & parser, bool check) { diff --git a/src/phreeqcpp/cxxKinetics.h b/src/phreeqcpp/cxxKinetics.h index 30555b98..9800e372 100644 --- a/src/phreeqcpp/cxxKinetics.h +++ b/src/phreeqcpp/cxxKinetics.h @@ -25,11 +25,6 @@ class cxxKinetics:public cxxNumKeyword void dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out=NULL) const; - void dump_essential_names(std::vector &e_names) const; - - void get_essential_values(std::vector &e_values) const; - void set_essential_values(std::vector::iterator &it); - void read_raw(CParser & parser, bool check = true); std::vector < LDBLE > &Get_steps(void) {return steps;} diff --git a/src/phreeqcpp/mainsubs.cpp b/src/phreeqcpp/mainsubs.cpp index baa20e71..fccfdf23 100644 --- a/src/phreeqcpp/mainsubs.cpp +++ b/src/phreeqcpp/mainsubs.cpp @@ -341,6 +341,107 @@ set_use(void) */ return (OK); } + +int Phreeqc::initial_solutions_poet(int sol_id) +/* ---------------------------------------------------------------------- */ +{ + /* + * Go through list of solutions, make initial solution calculations + * for any marked "new". + */ + int converge, converge1; + int last, n_user, print1; + char token[2 * MAX_LENGTH]; + + // state = INITIAL_SOLUTION; + set_use(); + print1 = TRUE; + dl_type_x = cxxSurface::NO_DL; + // std::map::iterator it = Rxn_solution_map.begin(); + // for ( ; it != Rxn_solution_map.end(); it++) + //{ + // for (size_t nn = 0; nn < Rxn_new_solution.size(); nn++) + cxxSolution &solution_ref = Rxn_solution_map.find(sol_id)->second; + initial_solution_isotopes = FALSE; + // if (solution_ref.Get_new_def()) { + use.Set_solution_ptr(&solution_ref); + LDBLE d0 = solution_ref.Get_density(); + // LDBLE d1 = 0; + bool diag = (diagonal_scale == TRUE) ? true : false; + int count_iterations = 0; + // std::string input_units = solution_ref.Get_initial_data()->Get_units(); + // cxxISolution *initial_data_ptr = solution_ref.Get_initial_data(); + density_iterations = 0; + // for (;;) { + prep(); + k_temp(solution_ref.Get_tc(), solution_ref.Get_patm()); + set(TRUE); + always_full_pitzer = FALSE; + + diagonal_scale = TRUE; + converge = model(); + if (converge == ERROR /*&& diagonal_scale == FALSE*/) { + diagonal_scale = TRUE; + always_full_pitzer = TRUE; + set(TRUE); + converge = model(); + } + // density_iterations++; + // if (solution_ref.Get_initial_data()->Get_calc_density()) { + // solution_ref.Set_density(calc_dens()); + // if (!equal(d0, solution_ref.Get_density(), 1e-8)) { + // initial_data_ptr->Set_units(input_units); + // d0 = solution_ref.Get_density(); + // if (count_iterations++ < 20) { + // diag = (diagonal_scale == TRUE) ? true : false; + // continue; + // } else { + // error_msg( + // sformatf("%s %d.", + // "Density calculation failed for initial solution ", + // solution_ref.Get_n_user()), + // STOP); + // } + // } + // } + // break; + // } + diagonal_scale = (diag) ? TRUE : FALSE; + converge1 = check_residuals(); + sum_species(); + viscosity(); + add_isotopes(solution_ref); + punch_all(); + // print_all(); + density_iterations = 0; + /* free_model_allocs(); */ + // remove pr_in + for (size_t i = 0; i < count_unknowns; i++) { + if (x[i]->type == SOLUTION_PHASE_BOUNDARY) + x[i]->phase->pr_in = false; + } + + if (converge == ERROR || converge1 == ERROR) { + error_msg(sformatf("%s %d.", + "Model failed to converge for initial solution ", + solution_ref.Get_n_user()), + STOP); + } + n_user = solution_ref.Get_n_user(); + last = solution_ref.Get_n_user_end(); + /* copy isotope data */ + if (solution_ref.Get_isotopes().size() > 0) { + isotopes_x = solution_ref.Get_isotopes(); + } else { + isotopes_x.clear(); + } + xsolution_save(n_user); + Utilities::Rxn_copies(Rxn_solution_map, n_user, last); + // } + initial_solution_isotopes = FALSE; + return (OK); +} + /* ---------------------------------------------------------------------- */ int Phreeqc:: initial_solutions(int print) diff --git a/src/phreeqcpp/model.cpp b/src/phreeqcpp/model.cpp index 60e077fb..103bcf4a 100644 --- a/src/phreeqcpp/model.cpp +++ b/src/phreeqcpp/model.cpp @@ -5003,16 +5003,16 @@ surface_model(void) */ debug_diffuse_layer_save = debug_diffuse_layer; debug_model_save = debug_model; - if (last_model.force_prep) - { + // if (last_model.force_prep) + // { + // same_model = FALSE; + // } + // else + // { + // same_model = check_same_model(); + // } same_model = FALSE; - } - else - { - same_model = check_same_model(); - } - if (dl_type_x != cxxSurface::NO_DL && same_model == FALSE) - { + if (dl_type_x != cxxSurface::NO_DL && same_model == FALSE) { s_diff_layer.clear(); for (i = 0; i < (int)s.size(); i++) { diff --git a/src/phreeqcpp/prep.cpp b/src/phreeqcpp/prep.cpp index a8349d3e..fade6839 100644 --- a/src/phreeqcpp/prep.cpp +++ b/src/phreeqcpp/prep.cpp @@ -36,21 +36,11 @@ prep(void) */ cxxSolution *solution_ptr; - if (state >= REACTION) - { - same_model = check_same_model(); - } - else - { - same_model = FALSE; - last_model.force_prep = true; - } - /*same_model = FALSE; */ -/* - * Initialize s, master, and unknown pointers - */ - solution_ptr = use.Get_solution_ptr(); - if (solution_ptr == NULL) + // UP: we force to reset the model + same_model = FALSE; + last_model.force_prep = true; + solution_ptr = use.Get_solution_ptr(); + if (solution_ptr == NULL) { error_msg("Solution needed for calculation not found, stopping.", STOP); @@ -87,12 +77,17 @@ prep(void) * Allocate space for array */ my_array.resize((max_unknowns + 1) * max_unknowns); - delta.resize(max_unknowns); - residual.resize(max_unknowns); - for (int j = 0; j < max_unknowns; j++) - { - residual[j] = 0; - } + for (double &val : my_array) { + val = 0; + } + delta.resize(max_unknowns); + for (double &val : delta) { + val = 0; + } + residual.resize(max_unknowns); + for (int j = 0; j < max_unknowns; j++) { + residual[j] = 0; + } /* * Build lists to fill Jacobian array and species list