[CL] Set aqueous solution.

ChemistryLib/PhreeqcIO.cpp

@@ -20,6 +20,7 @@
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTreeUtil.h"
 #include "MathLib/LinAlg/Eigen/EigenVector.h"
+#include "MathLib/LinAlg/LinAlg.h"
 #include "MeshLib/Mesh.h"
 #include "PhreeqcIOData/AqueousSolution.h"
 #include "PhreeqcIOData/ChemicalSystem.h"
@@ -116,8 +117,7 @@ void PhreeqcIO::initialize()
 void PhreeqcIO::executeInitialCalculation(
     std::vector<GlobalVector> const& interpolated_process_solutions)
 {
-    setAqueousSolutionsOrUpdateProcessSolutions(
-        process_solutions, Status::SettingAqueousSolutions);
+    setAqueousSolution(interpolated_process_solutions);
 
     writeInputsToFile();
 
@@ -132,8 +132,7 @@ void PhreeqcIO::executeInitialCalculation(
 void PhreeqcIO::doWaterChemistryCalculation(
     std::vector<GlobalVector*>& process_solutions, double const dt)
 {
-    setAqueousSolutionsOrUpdateProcessSolutions(
-        process_solutions, Status::SettingAqueousSolutions);
+    setAqueousSolution(interpolated_process_solutions);
 
     setAqueousSolutionsPrevFromDumpFile();
 
@@ -147,66 +146,23 @@ void PhreeqcIO::doWaterChemistryCalculation(
         process_solutions, Status::UpdatingProcessSolutions);
 }
 
-void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
-    std::vector<GlobalVector*> const& process_solutions, Status const status)
+void PhreeqcIO::setAqueousSolution(
+    std::vector<GlobalVector> const& interpolated_process_solutions)
 {
-    std::size_t const num_chemical_systems = _mesh.getNumberOfBaseNodes();
-
-    auto const chemical_system_map =
-        *_mesh.getProperties().template getPropertyVector<std::size_t>(
-            "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
-
-    // Loop over chemical systems
     auto& aqueous_solution = _chemical_system->aqueous_solution;
+
+    // components
     auto& components = aqueous_solution->components;
-    for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
+    for (unsigned component_id = 0; component_id < components.size();
+         ++component_id)
     {
-        auto const global_id = chemical_system_map[local_id];
-        // Loop over transport process id map to retrieve component
-        // concentrations from process solutions or to update process solutions
-        // after chemical calculation by Phreeqc
-        for (unsigned component_id = 0; component_id < components.size();
-             ++component_id)
-        {
-            auto& component = components[component_id];
-            auto& transport_process_solution =
-                process_solutions[component_id + 1];
-            switch (status)
-            {
-                case Status::SettingAqueousSolutions:
-                    // Set component concentrations.
-                    component.amount->set(
-                        local_id, (*transport_process_solution)[global_id]);
-                    break;
-                case Status::UpdatingProcessSolutions:
-                    // Update solutions of component transport processes.
-                    transport_process_solution->set(
-                        global_id, (*component.amount)[local_id]);
-                    break;
-            }
-        }
-
-        switch (status)
-        {
-            case Status::SettingAqueousSolutions:
-            {
-                // Set pH value by hydrogen concentration.
-                aqueous_solution->pH->set(
-                    local_id,
-                    -std::log10((*process_solutions.back())[global_id]));
-                break;
-            }
-            case Status::UpdatingProcessSolutions:
-            {
-                // Update hydrogen concentration by pH value.
-                auto hydrogen_concentration =
-                    std::pow(10, -(*aqueous_solution->pH)[local_id]);
-                process_solutions.back()->set(global_id,
-                                              hydrogen_concentration);
-                break;
-            }
-        }
+        MathLib::LinAlg::copy(interpolated_process_solutions[component_id + 1],
+                              *components[component_id].amount);
     }
+
+    // pH
+    MathLib::LinAlg::copy(interpolated_process_solutions.back(),
+                          *aqueous_solution->pH);
 }
 
 void PhreeqcIO::setAqueousSolutionsPrevFromDumpFile()

ChemistryLib/PhreeqcIO.h

@@ -31,12 +31,6 @@ struct SurfaceSite;
 struct Dump;
 struct UserPunch;
 
-enum class Status
-{
-    SettingAqueousSolutions,
-    UpdatingProcessSolutions
-};
-
 class PhreeqcIO final : public ChemicalSolverInterface
 {
 public:
@@ -60,9 +54,8 @@ public:
         std::vector<GlobalVector*>& process_solutions,
         double const dt) override;
 
-    void setAqueousSolutionsOrUpdateProcessSolutions(
-        std::vector<GlobalVector*> const& process_solutions,
-        Status const status);
+    void setAqueousSolution(
+        std::vector<GlobalVector> const& interpolated_process_solutions);
 
     void writeInputsToFile(double const dt = 0);