diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 17b1464d2484130833d0fb16a41d4a365e676208..f1d691ca3f73dc55a9cba4b39f6079920a493a4c 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -33,9 +33,15 @@ public:
return {};
}
+ virtual void computeSecondaryVariable(
+ std::size_t const /*ele_id*/,
+ std::vector<GlobalIndexType> const& /*chemical_system_indices*/)
+ {
+ }
+
virtual ~ChemicalSolverInterface() = default;
public:
- std::vector<std::vector<GlobalIndexType>> chemical_system_index_map;
+ std::vector<GlobalIndexType> chemical_system_index_map;
};
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 4f27d58b915e7507b10cfa277061f3091f26bf1d..3094cc50ec159616c328b4742801996742a5f7b8 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -104,11 +104,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name,
void PhreeqcIO::initialize()
{
- _num_chemical_systems = std::accumulate(
- begin(chemical_system_index_map), end(chemical_system_index_map), 0,
- [](auto result, auto const& chemical_system_index) {
- return result + chemical_system_index.size();
- });
+ _num_chemical_systems = chemical_system_index_map.size();
_chemical_system->initialize(_num_chemical_systems);
@@ -248,83 +244,73 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
os << reaction_rates << "\n";
}
- for (std::size_t mesh_item_id = 0;
- mesh_item_id < phreeqc_io.chemical_system_index_map.size();
- ++mesh_item_id)
+ for (std::size_t chemical_system_id = 0;
+ chemical_system_id < phreeqc_io._num_chemical_systems;
+ ++chemical_system_id)
{
- auto const& local_chemical_system =
- phreeqc_io.chemical_system_index_map[mesh_item_id];
- for (std::size_t ip = 0; ip < local_chemical_system.size(); ++ip)
- {
- auto const& chemical_system_id = local_chemical_system[ip];
-
- os << "SOLUTION " << chemical_system_id + 1 << "\n";
- phreeqc_io._chemical_system->aqueous_solution->print(
- os, chemical_system_id);
+ os << "SOLUTION " << chemical_system_id + 1 << "\n";
+ phreeqc_io._chemical_system->aqueous_solution->print(
+ os, chemical_system_id);
- auto const& dump = phreeqc_io._dump;
- if (dump)
+ auto const& dump = phreeqc_io._dump;
+ if (dump)
+ {
+ auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
+ if (!aqueous_solutions_prev.empty())
{
- auto const& aqueous_solutions_prev =
- dump->aqueous_solutions_prev;
- if (!aqueous_solutions_prev.empty())
- {
- os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
- }
+ os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
}
+ }
- os << "USE solution none"
- << "\n";
- os << "END"
- << "\n\n";
+ os << "USE solution none"
+ << "\n";
+ os << "END"
+ << "\n\n";
- os << "USE solution " << chemical_system_id + 1 << "\n\n";
+ os << "USE solution " << chemical_system_id + 1 << "\n\n";
- auto const& equilibrium_reactants =
- phreeqc_io._chemical_system->equilibrium_reactants;
- if (!equilibrium_reactants.empty())
- {
- os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
- for (auto const& equilibrium_reactant : equilibrium_reactants)
- {
- equilibrium_reactant.print(os, chemical_system_id);
- }
- os << "\n";
- }
-
- auto const& kinetic_reactants =
- phreeqc_io._chemical_system->kinetic_reactants;
- if (!kinetic_reactants.empty())
+ auto const& equilibrium_reactants =
+ phreeqc_io._chemical_system->equilibrium_reactants;
+ if (!equilibrium_reactants.empty())
+ {
+ os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
+ for (auto const& equilibrium_reactant : equilibrium_reactants)
{
- os << "KINETICS " << chemical_system_id + 1 << "\n";
- for (auto const& kinetic_reactant : kinetic_reactants)
- {
- kinetic_reactant.print(os, chemical_system_id);
- }
- os << "-steps " << phreeqc_io._dt << "\n"
- << "\n";
+ equilibrium_reactant.print(os, chemical_system_id);
}
+ os << "\n";
+ }
- auto const& surface = phreeqc_io._surface;
- if (!surface.empty())
+ auto const& kinetic_reactants =
+ phreeqc_io._chemical_system->kinetic_reactants;
+ if (!kinetic_reactants.empty())
+ {
+ os << "KINETICS " << chemical_system_id + 1 << "\n";
+ for (auto const& kinetic_reactant : kinetic_reactants)
{
- os << "SURFACE " << chemical_system_id + 1 << "\n";
- std::size_t aqueous_solution_id =
- dump->aqueous_solutions_prev.empty()
- ? chemical_system_id + 1
- : phreeqc_io._num_chemical_systems +
- chemical_system_id + 1;
- os << "-equilibrate with solution " << aqueous_solution_id
- << "\n";
- os << "-sites_units DENSITY"
- << "\n";
- os << surface << "\n";
- os << "SAVE solution " << chemical_system_id + 1 << "\n";
+ kinetic_reactant.print(os, chemical_system_id);
}
+ os << "-steps " << phreeqc_io._dt << "\n"
+ << "\n";
+ }
- os << "END"
- << "\n\n";
+ auto const& surface = phreeqc_io._surface;
+ if (!surface.empty())
+ {
+ os << "SURFACE " << chemical_system_id + 1 << "\n";
+ std::size_t aqueous_solution_id =
+ dump->aqueous_solutions_prev.empty()
+ ? chemical_system_id + 1
+ : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
+ os << "-equilibrate with solution " << aqueous_solution_id << "\n";
+ os << "-sites_units DENSITY"
+ << "\n";
+ os << surface << "\n";
+ os << "SAVE solution " << chemical_system_id + 1 << "\n";
}
+
+ os << "END"
+ << "\n\n";
}
auto const& dump = phreeqc_io._dump;
@@ -387,186 +373,145 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
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)
- {
- std::fill(kinetic_reactant.amount_avg->begin(),
- kinetic_reactant.amount_avg->end(), 0.);
- }
- for (std::size_t mesh_item_id = 0;
- mesh_item_id < phreeqc_io.chemical_system_index_map.size();
- ++mesh_item_id)
+ for (std::size_t chemical_system_id = 0;
+ chemical_system_id < phreeqc_io._num_chemical_systems;
+ ++chemical_system_id)
{
- auto const& local_chemical_system =
- phreeqc_io.chemical_system_index_map[mesh_item_id];
- auto const num_local_chemical_system = local_chemical_system.size();
- for (std::size_t ip = 0; ip < num_local_chemical_system; ++ip)
+ // Skip equilibrium calculation result of initial solution
+ for (int i = 0; i < num_skipped_lines; ++i)
{
- auto const& chemical_system_id = local_chemical_system[ip];
-
- // Skip equilibrium calculation result of initial solution
- for (int i = 0; i < num_skipped_lines; ++i)
- {
- in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
- }
+ in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+ }
- // Get calculation result of the solution after the reaction
- if (!std::getline(in, line))
- {
- OGS_FATAL(
- "Error when reading calculation result of Solution {:d} "
- "after "
- "the reaction.",
- chemical_system_id);
- }
+ // Get calculation result of the solution after the reaction
+ if (!std::getline(in, line))
+ {
+ OGS_FATAL(
+ "Error when reading calculation result of Solution {:d} "
+ "after "
+ "the reaction.",
+ chemical_system_id);
+ }
- std::vector<double> accepted_items;
- std::vector<std::string> items;
- boost::trim_if(line, boost::is_any_of("\t "));
- boost::algorithm::split(items, line, boost::is_any_of("\t "),
- boost::token_compress_on);
- for (int item_id = 0; item_id < static_cast<int>(items.size());
- ++item_id)
+ std::vector<double> accepted_items;
+ std::vector<std::string> items;
+ boost::trim_if(line, boost::is_any_of("\t "));
+ boost::algorithm::split(items, line, boost::is_any_of("\t "),
+ boost::token_compress_on);
+ for (int item_id = 0; item_id < static_cast<int>(items.size());
+ ++item_id)
+ {
+ if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
+ item_id) == dropped_item_ids.end())
{
- if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
- item_id) == dropped_item_ids.end())
+ double value;
+ try
{
- double value;
- try
- {
- value = std::stod(items[item_id]);
- }
- catch (const std::invalid_argument& e)
- {
- OGS_FATAL(
- "Invalid argument. Could not convert string '{:s}' "
- "to double for chemical system {:d}, column {:d}. "
- "Exception '{:s}' was thrown.",
- items[item_id], chemical_system_id + 1, item_id,
- e.what());
- }
- catch (const std::out_of_range& e)
- {
- OGS_FATAL(
- "Out of range error. Could not convert string "
- "'{:s}' to double for chemical system {:d}, column "
- "{:d}. Exception '{:s}' was thrown.",
- items[item_id], chemical_system_id + 1, item_id,
- e.what());
- }
- accepted_items.push_back(value);
+ value = std::stod(items[item_id]);
}
- }
- assert(accepted_items.size() == output.accepted_items.size());
-
- auto& aqueous_solution =
- phreeqc_io._chemical_system->aqueous_solution;
- auto& components = aqueous_solution->components;
- auto& user_punch = phreeqc_io._user_punch;
- for (int item_id = 0;
- item_id < static_cast<int>(accepted_items.size());
- ++item_id)
- {
- auto const& accepted_item = output.accepted_items[item_id];
- auto const& item_name = accepted_item.name;
-
- auto compare_by_name = [&item_name](auto const& item) {
- return item.name == item_name;
- };
-
- switch (accepted_item.item_type)
+ catch (const std::invalid_argument& e)
{
- case ItemType::pH:
- {
- // Update pH value
- aqueous_solution->pH->set(
- chemical_system_id,
- std::pow(10, -accepted_items[item_id]));
- break;
- }
- case ItemType::pe:
- {
- // Update pe value
- (*aqueous_solution->pe)[chemical_system_id] =
- accepted_items[item_id];
- break;
- }
- case ItemType::Component:
- {
- // Update component concentrations
- auto& component = BaseLib::findElementOrError(
- components.begin(), components.end(),
- compare_by_name,
- "Could not find component '" + item_name + "'.");
- component.amount->set(chemical_system_id,
- accepted_items[item_id]);
- break;
- }
- case ItemType::EquilibriumReactant:
- {
- // Update amounts of equilibrium reactant
- auto& equilibrium_reactant =
- BaseLib::findElementOrError(
- equilibrium_reactants.begin(),
- equilibrium_reactants.end(), compare_by_name,
- "Could not find equilibrium reactant '" +
- 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:
- {
- // Update amounts of kinetic reactants
- auto& kinetic_reactant = BaseLib::findElementOrError(
- kinetic_reactants.begin(), kinetic_reactants.end(),
- compare_by_name,
- "Could not find kinetic reactant '" + item_name +
- "'.");
- (*kinetic_reactant.amount)[chemical_system_id] =
- accepted_items[item_id];
- (*kinetic_reactant.amount_avg)[mesh_item_id] +=
- accepted_items[item_id];
- break;
- }
- case ItemType::SecondaryVariable:
- {
- assert(user_punch);
- auto& secondary_variables =
- user_punch->secondary_variables;
- // Update values of secondary variables
- auto& secondary_variable = BaseLib::findElementOrError(
- secondary_variables.begin(),
- secondary_variables.end(), compare_by_name,
- "Could not find secondary variable '" + item_name +
- "'.");
- (*secondary_variable.value)[chemical_system_id] =
- accepted_items[item_id];
- break;
- }
+ OGS_FATAL(
+ "Invalid argument. Could not convert string '{:s}' "
+ "to double for chemical system {:d}, column {:d}. "
+ "Exception '{:s}' was thrown.",
+ items[item_id], chemical_system_id + 1, item_id,
+ e.what());
}
+ catch (const std::out_of_range& e)
+ {
+ OGS_FATAL(
+ "Out of range error. Could not convert string "
+ "'{:s}' to double for chemical system {:d}, column "
+ "{:d}. Exception '{:s}' was thrown.",
+ items[item_id], chemical_system_id + 1, item_id,
+ e.what());
+ }
+ accepted_items.push_back(value);
}
}
+ assert(accepted_items.size() == output.accepted_items.size());
- for (auto& equilibrium_reactant : equilibrium_reactants)
+ auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
+ auto& components = aqueous_solution->components;
+ auto& user_punch = phreeqc_io._user_punch;
+ for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
+ ++item_id)
{
- (*equilibrium_reactant.amount_avg)[mesh_item_id] /=
- num_local_chemical_system;
- }
+ auto const& accepted_item = output.accepted_items[item_id];
+ auto const& item_name = accepted_item.name;
- for (auto& kinetic_reactant : kinetic_reactants)
- {
- (*kinetic_reactant.amount_avg)[mesh_item_id] /=
- num_local_chemical_system;
+ auto compare_by_name = [&item_name](auto const& item) {
+ return item.name == item_name;
+ };
+
+ switch (accepted_item.item_type)
+ {
+ case ItemType::pH:
+ {
+ // Update pH value
+ aqueous_solution->pH->set(
+ chemical_system_id,
+ std::pow(10, -accepted_items[item_id]));
+ break;
+ }
+ case ItemType::pe:
+ {
+ // Update pe value
+ (*aqueous_solution->pe)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ case ItemType::Component:
+ {
+ // Update component concentrations
+ auto& component = BaseLib::findElementOrError(
+ components.begin(), components.end(), compare_by_name,
+ "Could not find component '" + item_name + "'.");
+ component.amount->set(chemical_system_id,
+ accepted_items[item_id]);
+ break;
+ }
+ case ItemType::EquilibriumReactant:
+ {
+ // Update amounts of equilibrium reactant
+ auto& equilibrium_reactant = BaseLib::findElementOrError(
+ equilibrium_reactants.begin(),
+ equilibrium_reactants.end(), compare_by_name,
+ "Could not find equilibrium reactant '" + item_name +
+ "'.");
+ (*equilibrium_reactant.amount)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ case ItemType::KineticReactant:
+ {
+ // Update amounts of kinetic reactants
+ auto& kinetic_reactant = BaseLib::findElementOrError(
+ kinetic_reactants.begin(), kinetic_reactants.end(),
+ compare_by_name,
+ "Could not find kinetic reactant '" + item_name + "'.");
+ (*kinetic_reactant.amount)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ case ItemType::SecondaryVariable:
+ {
+ assert(user_punch);
+ auto& secondary_variables = user_punch->secondary_variables;
+ // Update values of secondary variables
+ auto& secondary_variable = BaseLib::findElementOrError(
+ secondary_variables.begin(), secondary_variables.end(),
+ compare_by_name,
+ "Could not find secondary variable '" + item_name +
+ "'.");
+ (*secondary_variable.value)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ }
}
}
@@ -585,5 +530,35 @@ std::vector<std::string> const PhreeqcIO::getComponentList() const
return component_names;
}
+
+template <typename Reactant>
+static double averageReactantAmount(
+ Reactant const& reactant,
+ std::vector<GlobalIndexType> const& chemical_system_indices)
+{
+ double const sum = std::accumulate(
+ chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
+ [&](double const s, GlobalIndexType const id) {
+ return s + (*reactant.amount)[id];
+ });
+ return sum / chemical_system_indices.size();
+}
+
+void PhreeqcIO::computeSecondaryVariable(
+ std::size_t const ele_id,
+ std::vector<GlobalIndexType> const& chemical_system_indices)
+{
+ for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+ {
+ (*kinetic_reactant.amount_avg)[ele_id] =
+ averageReactantAmount(kinetic_reactant, chemical_system_indices);
+ }
+
+ for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+ {
+ (*equilibrium_reactant.amount_avg)[ele_id] = averageReactantAmount(
+ equilibrium_reactant, chemical_system_indices);
+ }
+}
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index ab15c9d061074e7534d34e75f9d8fd9729abcdae..fb7b20ef63dccc11fe599dde2160ea169b613266 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;
@@ -92,7 +96,7 @@ private:
Knobs const _knobs;
double _dt = std::numeric_limits<double>::quiet_NaN();
const int phreeqc_instance_id = 0;
- int _num_chemical_systems = -1;
+ std::size_t _num_chemical_systems = -1;
};
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ProcessLib/ComponentTransport/ChemicalProcessData.h b/ProcessLib/ComponentTransport/ChemicalProcessData.h
deleted file mode 100644
index 8b805302324f0bbbc2332bfa878554f17c1d8bd8..0000000000000000000000000000000000000000
--- a/ProcessLib/ComponentTransport/ChemicalProcessData.h
+++ /dev/null
@@ -1,32 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <vector>
-
-#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-struct ChemicalProcessData
-{
- ChemicalProcessData(
- std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map_)
- : chemical_system_index_map(chemical_system_index_map_)
- {
- }
-
- std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map;
-};
-} // namespace ComponentTransport
-} // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 42524cd0a607d8ab5a1ecb58aab87b5dcb820786..c0f79544127afc9c7c7a5fc75a7a4be71590b50e 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -48,6 +48,8 @@ struct IntegrationPointData final
GlobalDimNodalMatrixType const dNdx;
double const integration_weight;
+ GlobalIndexType chemical_system_id = 0;
+
double porosity = std::numeric_limits<double>::quiet_NaN();
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
@@ -66,6 +68,8 @@ public:
_coupled_solutions = coupling_term;
}
+ virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
+
virtual std::vector<double> const& getIntPtDarcyVelocity(
const double t,
std::vector<GlobalVector*> const& x,
@@ -159,23 +163,23 @@ public:
medium->property(MaterialPropertyLib::PropertyType::porosity)
.template initialValue<double>(pos, 0.);
}
+ }
- if (_process_data.chemical_process_data)
- {
- // chemical system index map
- auto& chemical_system_index_map =
- _process_data.chemical_process_data->chemical_system_index_map;
-
- GlobalIndexType const start_value =
- chemical_system_index_map.empty()
- ? 0
- : chemical_system_index_map.back().back() + 1;
-
- std::vector<GlobalIndexType> indices(
- _integration_method.getNumberOfPoints());
- std::iota(indices.begin(), indices.end(), start_value);
+ void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
+ {
+ assert(_process_data.chemical_solver_interface);
+ // chemical system index map
+ auto& chemical_system_index_map =
+ _process_data.chemical_solver_interface->chemical_system_index_map;
- chemical_system_index_map.push_back(std::move(indices));
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data[ip].chemical_system_id = chemical_system_index_map.empty()
+ ? 0
+ : chemical_system_index_map.back() + 1;
+ chemical_system_index_map.push_back(_ip_data[ip].chemical_system_id);
}
}
@@ -945,14 +949,19 @@ public:
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
- assert(_process_data.chemical_process_data);
+ assert(_process_data.chemical_solver_interface);
assert(int_pt_x.size() == 1);
cache.clear();
- auto const ele_id = _element.getID();
- auto const& indices = _process_data.chemical_process_data
- ->chemical_system_index_map[ele_id];
- cache = int_pt_x[0]->get(indices);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip(0); ip < n_integration_points; ++ip)
+ {
+ auto const& chemical_system_id = _ip_data[ip].chemical_system_id;
+ auto const c_int_pt = int_pt_x[0]->get(chemical_system_id);
+ cache.push_back(c_int_pt);
+ }
return cache;
}
@@ -976,11 +985,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*/,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index b0f217f60f59791863e8e1b589c5ccac00550aeb..f6c880a3537eaf39f10afaf09a7dbfd8c81c18e4 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -86,6 +86,10 @@ void ComponentTransportProcess::initializeConcreteProcess(
if (_chemical_solver_interface)
{
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &ComponentTransportLocalAssemblerInterface::setChemicalSystemID,
+ _local_assemblers, pv.getActiveElementIDs());
+
_chemical_solver_interface->initialize();
}
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
index 25f9c386e5b4181d2bfce46208c2568eec61619a..bc95bca0dd1d420f526bfd4166109ba9313efad3 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
@@ -12,7 +12,7 @@
#include <memory>
-#include "ChemicalProcessData.h"
+#include "ChemistryLib/ChemicalSolverInterface.h"
#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
@@ -38,7 +38,7 @@ struct ComponentTransportProcessData
Eigen::VectorXd const specific_body_force;
bool const has_gravity;
bool const non_advective_form;
- std::unique_ptr<ChemicalProcessData> chemical_process_data;
+ ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface;
const int hydraulic_process_id;
// TODO (renchao-lu): This variable is used in the calculation of the
diff --git a/ProcessLib/ComponentTransport/CreateChemicalProcessData.cpp b/ProcessLib/ComponentTransport/CreateChemicalProcessData.cpp
deleted file mode 100644
index d60dc12ca414a235ce078c613ccd7538830b1ff4..0000000000000000000000000000000000000000
--- a/ProcessLib/ComponentTransport/CreateChemicalProcessData.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "CreateChemicalProcessData.h"
-#include "ChemicalProcessData.h"
-
-#include "ChemistryLib/ChemicalSolverInterface.h"
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-std::unique_ptr<ChemicalProcessData> createChemicalProcessData(
- ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface)
-{
- if (!chemical_solver_interface)
- {
- return nullptr;
- }
-
- return std::make_unique<ChemicalProcessData>(
- chemical_solver_interface->chemical_system_index_map);
-}
-} // namespace ComponentTransport
-} // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/CreateChemicalProcessData.h b/ProcessLib/ComponentTransport/CreateChemicalProcessData.h
deleted file mode 100644
index bd8b3e491250dbdca6bad0122ba4d779989035d7..0000000000000000000000000000000000000000
--- a/ProcessLib/ComponentTransport/CreateChemicalProcessData.h
+++ /dev/null
@@ -1,29 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <memory>
-
-namespace ChemistryLib
-{
-class ChemicalSolverInterface;
-}
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-struct ChemicalProcessData;
-
-std::unique_ptr<ChemicalProcessData> createChemicalProcessData(
- ChemistryLib::ChemicalSolverInterface* chemical_solver_interface);
-} // namespace ComponentTransport
-} // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index d04aa117f4a9ee633592ccd62182023379deb288..c025545505cb2e37cdf74d24cfd02b66a6654c47 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -9,7 +9,6 @@
*/
#include "CreateComponentTransportProcess.h"
-#include "CreateChemicalProcessData.h"
#include "ChemistryLib/ChemicalSolverInterface.h"
#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
@@ -225,14 +224,11 @@ std::unique_ptr<Process> createComponentTransportProcess(
checkMPLProperties(mesh, *media_map);
DBUG("Media properties verified.");
- auto chemical_process_data =
- createChemicalProcessData(chemical_solver_interface.get());
-
ComponentTransportProcessData process_data{std::move(media_map),
specific_body_force,
has_gravity,
non_advective_form,
- std::move(chemical_process_data),
+ chemical_solver_interface.get(),
hydraulic_process_id,
first_transport_process_id};