From ef39384a3b63cca39d8ccab036e817f8134dd426 Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Wed, 23 Dec 2020 12:38:24 +0100
Subject: [PATCH] [PL] Calculate averaged solid mass for output.
---
ChemistryLib/ChemicalSolverInterface.h | 6 ++++
ChemistryLib/PhreeqcIO.cpp | 30 +++++++++++++++++++
ChemistryLib/PhreeqcIO.h | 4 +++
.../ComponentTransportFEM.h | 15 +++++++++-
4 files changed, 54 insertions(+), 1 deletion(-)
diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 21d53e9fd14..f1d691ca3f7 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/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:
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 45af3e26944..8149ec0a32e 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/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
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 6df10d8efbe..fb7b20ef63d 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/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;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 957e5c8bbc8..27bb0d6c974 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/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*/,
--
GitLab