[CL] Output equilibrium reactants on the element level.

ChemistryLib/PhreeqcIO.cpp

@@ -373,6 +373,14 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
     auto const& surface = phreeqc_io._surface;
     int const num_skipped_lines = surface.empty() ? 1 : 2;
 
+    auto& equilibrium_reactants =
+        phreeqc_io._chemical_system->equilibrium_reactants;
+    for (auto& equilibrium_reactant : equilibrium_reactants)
+    {
+        std::fill(equilibrium_reactant.amount_avg->begin(),
+                  equilibrium_reactant.amount_avg->end(), 0.);
+    }
+
     auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
     for (auto& kinetic_reactant : kinetic_reactants)
     {
@@ -449,8 +457,6 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
             auto& aqueous_solution =
                 phreeqc_io._chemical_system->aqueous_solution;
             auto& components = aqueous_solution->components;
-            auto& equilibrium_reactants =
-                phreeqc_io._chemical_system->equilibrium_reactants;
             auto& user_punch = phreeqc_io._user_punch;
             for (int item_id = 0;
                  item_id < static_cast<int>(accepted_items.size());
@@ -502,6 +508,8 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                                     item_name + "'.");
                         (*equilibrium_reactant.amount)[chemical_system_id] =
                             accepted_items[item_id];
+                        (*equilibrium_reactant.amount_avg)[mesh_item_id] +=
+                            accepted_items[item_id];
                         break;
                     }
                     case ItemType::KineticReactant:
@@ -537,6 +545,12 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
             }
         }
 
+        for (auto& equilibrium_reactant : equilibrium_reactants)
+        {
+            (*equilibrium_reactant.amount_avg)[mesh_item_id] /=
+                num_local_chemical_system;
+        }
+
         for (auto& kinetic_reactant : kinetic_reactants)
         {
             (*kinetic_reactant.amount_avg)[mesh_item_id] /=

ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp

@@ -50,8 +50,15 @@ std::vector<EquilibriumReactant> createEquilibriumReactants(
         auto amount = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
-        equilibrium_reactants.emplace_back(
-            std::move(name), amount, initial_amount, saturation_index);
+        auto mesh_prop_amount = MeshLib::getOrCreateMeshProperty<double>(
+            mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
+        mesh_prop_amount->resize(mesh.getNumberOfElements());
+
+        equilibrium_reactants.emplace_back(std::move(name),
+                                           amount,
+                                           mesh_prop_amount,
+                                           initial_amount,
+                                           saturation_index);
     }
 
     return equilibrium_reactants;

ChemistryLib/PhreeqcIOData/EquilibriumReactant.h

@@ -30,10 +30,12 @@ struct EquilibriumReactant
 {
     EquilibriumReactant(std::string name_,
                         MeshLib::PropertyVector<double>* amount_,
+                        MeshLib::PropertyVector<double>* amount_avg_,
                         double const initial_amount_,
                         double saturation_index_)
         : name(std::move(name_)),
           amount(amount_),
+          amount_avg(amount_avg_),
           initial_amount(initial_amount_),
           saturation_index(saturation_index_)
     {
@@ -43,6 +45,7 @@ struct EquilibriumReactant
 
     std::string const name;
     MeshLib::PropertyVector<double>* amount;
+    MeshLib::PropertyVector<double>* amount_avg;
     double const initial_amount;
     double const saturation_index;
     static const ItemType item_type = ItemType::EquilibriumReactant;