diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index a1a6eb1f66ee705e067fb229b03f64b24ad586d2..79443449e6540fafb40ef956058c3a44cc1c47a4 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -30,10 +30,16 @@ public:
virtual void initialize() {}
virtual void initializeChemicalSystemConcrete(
+ std::vector<double> const& /*concentrations*/,
GlobalIndexType const& /*chemical_system_id*/,
MaterialPropertyLib::Medium const* /*medium*/,
- ParameterLib::SpatialPosition const& /*pos*/,
- double const /*t*/)
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/)
+ {
+ }
+
+ virtual void setChemicalSystemConcrete(
+ std::vector<double> const& /*concentrations*/,
+ GlobalIndexType const& /*chemical_system_id*/)
{
}
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 1d77d3e1edd4e4539a866a40a27e4cf93709eef2..3ec32c7ee4f2f3006a6d7f372fca6996567b6f7b 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -50,6 +50,23 @@ std::ostream& operator<<(std::ostream& os,
return os;
}
+void setAqueousSolution(std::vector<double> const& concentrations,
+ GlobalIndexType const& chemical_system_id,
+ AqueousSolution& aqueous_solution)
+{
+ // components
+ auto& components = aqueous_solution.components;
+ for (unsigned component_id = 0; component_id < components.size();
+ ++component_id)
+ {
+ components[component_id].amount->set(chemical_system_id,
+ concentrations[component_id]);
+ }
+
+ // pH
+ aqueous_solution.pH->set(chemical_system_id, concentrations.back());
+}
+
template <typename Reactant>
void setReactantMolality(Reactant& reactant,
GlobalIndexType const& chemical_system_id,
@@ -145,10 +162,15 @@ void PhreeqcIO::initialize()
}
void PhreeqcIO::initializeChemicalSystemConcrete(
+ std::vector<double> const& concentrations,
GlobalIndexType const& chemical_system_id,
MaterialPropertyLib::Medium const* medium,
- ParameterLib::SpatialPosition const& pos, double const t)
+ ParameterLib::SpatialPosition const& pos,
+ double const t)
{
+ setAqueousSolution(concentrations, chemical_system_id,
+ *_chemical_system->aqueous_solution);
+
auto const& solid_phase = medium->phase("Solid");
for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
{
@@ -163,7 +185,13 @@ void PhreeqcIO::initializeChemicalSystemConcrete(
}
}
+void PhreeqcIO::setChemicalSystemConcrete(
+ std::vector<double> const& concentrations,
+ GlobalIndexType const& chemical_system_id)
{
+ setAqueousSolution(concentrations, chemical_system_id,
+ *_chemical_system->aqueous_solution);
+}
void PhreeqcIO::executeInitialCalculation()
{
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 44e6814f27d47fcd025ba87c70f5e53d9e724303..1e9d189ba94e0804376381c11f395e62eb6bbf2f 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -47,11 +47,15 @@ public:
void initialize() override;
void initializeChemicalSystemConcrete(
+ std::vector<double> const& concentrations,
GlobalIndexType const& chemical_system_id,
MaterialPropertyLib::Medium const* medium,
ParameterLib::SpatialPosition const& pos,
double const t) override;
+ void setChemicalSystemConcrete(
+ std::vector<double> const& concentrations,
+ GlobalIndexType const& chemical_system_id) override;
void executeInitialCalculation() override;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 86ffe60394d8679000cfe9ddd0fb999652b32b0f..b117cab62802e0289a9e896a488d30b738b3f4f5 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -97,6 +97,29 @@ public:
initializeChemicalSystemConcrete(local_x, t);
}
+ void setChemicalSystem(
+ std::size_t const mesh_item_id,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
+ std::vector<GlobalVector*> const& x, double const t, double const dt)
+ {
+ std::vector<double> local_x_vec;
+
+ auto const n_processes = x.size();
+ for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
+ {
+ auto const indices =
+ NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
+ assert(!indices.empty());
+ auto const local_solution = x[process_id]->get(indices);
+ local_x_vec.insert(std::end(local_x_vec),
+ std::begin(local_solution),
+ std::end(local_solution));
+ }
+ auto const local_x = MathLib::toVector(local_x_vec);
+
+ setChemicalSystemConcrete(local_x, t, dt);
+ }
+
virtual std::vector<double> const& getIntPtDarcyVelocity(
const double t,
std::vector<GlobalVector*> const& x,
@@ -113,6 +136,10 @@ private:
virtual void initializeChemicalSystemConcrete(
Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
+ virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
+ double const /*t*/,
+ double const /*dt*/) = 0;
+
protected:
CoupledSolutionsForStaggeredScheme* _coupled_solutions{nullptr};
};
@@ -217,7 +244,7 @@ public:
}
}
- void initializeChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
+ void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x,
double const t) override
{
assert(_process_data.chemical_solver_interface);
@@ -233,11 +260,62 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
auto const& chemical_system_id = ip_data.chemical_system_id;
+ auto const n_component = _transport_process_variables.size();
+ std::vector<double> C_int_pt(n_component);
+ for (unsigned component_id = 0; component_id < n_component;
+ ++component_id)
+ {
+ auto const concentration_index =
+ first_concentration_index +
+ component_id * concentration_size;
+ auto const local_C =
+ local_x.template segment<concentration_size>(
+ concentration_index);
+
+ NumLib::shapeFunctionInterpolate(local_C, N,
+ C_int_pt[component_id]);
+ }
+
_process_data.chemical_solver_interface
- ->initializeChemicalSystemConcrete(chemical_system_id, medium,
- pos, t);
+ ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
+ medium, pos, t);
+ }
+ }
+
+ void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
+ double const /*t*/, double /*dt*/) override
+ {
+ assert(_process_data.chemical_solver_interface);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& chemical_system_id = ip_data.chemical_system_id;
+
+ auto const n_component = _transport_process_variables.size();
+ std::vector<double> C_int_pt(n_component);
+ for (unsigned component_id = 0; component_id < n_component;
+ ++component_id)
+ {
+ auto const concentration_index =
+ first_concentration_index +
+ component_id * concentration_size;
+ auto const local_C =
+ local_x.template segment<concentration_size>(
+ concentration_index);
+
+ NumLib::shapeFunctionInterpolate(local_C, N,
+ C_int_pt[component_id]);
+ }
+
+ _process_data.chemical_solver_interface->setChemicalSystemConcrete(
+ C_int_pt, chemical_system_id);
}
}
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 5037f7a70b065d92115ae29441adde5100680370..8b59ddfecddb3de5a8ed42e23caecf9c68dcfa77 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -224,6 +224,18 @@ void ComponentTransportProcess::solveReactionEquation(
// process.
// TODO: move into a global loop to consider both mass balance over
// space and localized chemical equilibrium between solutes.
+ const int process_id = 0;
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+ dof_tables.reserve(x.size());
+ std::generate_n(std::back_inserter(dof_tables), x.size(),
+ [&]() { return _local_to_global_index_map.get(); });
+
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &ComponentTransportLocalAssemblerInterface::setChemicalSystem,
+ _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt);
+
BaseLib::RunTime time_phreeqc;
time_phreeqc.start();