[PL] Calculate averaged solid mass for output.

ChemistryLib/ChemicalSolverInterface.h

@@ -33,6 +33,12 @@ public:
         return {};
     }
 
+    virtual void computeSecondaryVariable(
+        std::size_t const /*ele_id*/,
+        std::vector<GlobalIndexType> const& /*chemical_system_indices*/)
+    {
+    }
+
     virtual ~ChemicalSolverInterface() = default;
 
 public:

ChemistryLib/PhreeqcIO.cpp

@@ -530,5 +530,35 @@ std::vector<std::string> const PhreeqcIO::getComponentList() const
 
     return component_names;
 }
+
+void PhreeqcIO::computeSecondaryVariable(
+    std::size_t const ele_id,
+    std::vector<GlobalIndexType> const& chemical_system_indices)
+{
+    auto const num_chemical_systems = chemical_system_indices.size();
+    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+    {
+        double amount_avg = std::accumulate(
+            chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
+            [&](double amount, auto const& id) {
+                return amount + (*kinetic_reactant.amount)[id];
+            });
+
+        amount_avg /= num_chemical_systems;
+        (*kinetic_reactant.amount_avg)[ele_id] = amount_avg;
+    }
+
+    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+    {
+        double amount_avg = std::accumulate(
+            chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
+            [&](double amount, auto const& id) {
+                return amount + (*equilibrium_reactant.amount)[id];
+            });
+
+        amount_avg /= num_chemical_systems;
+        (*equilibrium_reactant.amount_avg)[ele_id] = amount_avg;
+    }
+}
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib

ChemistryLib/PhreeqcIO.h

@@ -69,6 +69,10 @@ public:
 
     friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
 
+    void computeSecondaryVariable(
+        std::size_t const ele_id,
+        std::vector<GlobalIndexType> const& chemical_system_indices) override;
+
     std::vector<std::string> const getComponentList() const override;
 
     std::string const _phreeqc_input_file;

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -980,11 +980,24 @@ public:
         auto const ele_velocity_mat =
             MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
 
-        auto ele_id = _element.getID();
+        auto const ele_id = _element.getID();
         Eigen::Map<LocalVectorType>(
             &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
             GlobalDim) =
             ele_velocity_mat.rowwise().sum() / n_integration_points;
+
+        if (_process_data.chemical_solver_interface)
+        {
+            std::vector<GlobalIndexType> chemical_system_indices;
+            chemical_system_indices.reserve(n_integration_points);
+            std::transform(
+                _ip_data.begin(), _ip_data.end(),
+                std::back_inserter(chemical_system_indices),
+                [](auto const& ip_data) { return ip_data.chemical_system_id; });
+
+            _process_data.chemical_solver_interface->computeSecondaryVariable(
+                ele_id, chemical_system_indices);
+        }
     }
 
     void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,